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general-relativistic magnetohydrodynamics in arbitrary spacetimes and exploits 


adaptive mesh refinement techniques with an efficient block-based approach. 
Several spacetimes have already been implemented and tested. We demonstrate 
the validity of BHAC by means of various one-, two-, and three-dimensional test 
problems, as well as through a close comparison with the HARM3D code in the 
case of a torus accreting onto a black hole. The convergence of a turbulent 
accretion scenario is investigated with several diagnostics and we find accretion 
rates and horizon-penetrating fluxes to be convergent to within a few percent 
when the problem is run in three dimensions. Our analysis also involves the study 
of the corresponding thermal synchrotron emission, which is performed by means 
of a new general-relativistic radiative transfer code, BHOSS The resulting synthetic 
intensity maps of accretion onto black holes are found to be convergent with 
increasing resolution and are anticipated to play a crucial role in the 
interpretation of horizon-scale images resulting from upcoming radio observations 
of the source at the Galactic Center. 


1 Introduction 
Accreting black holes (BHs) are amongst the most powerful astrophysical objects 
in the Universe. A substantial fraction of the gravitational binding energy of the 
accreting gas is released within tens of gravitational radii from the BH, and this en- 
ergy supplies the power for a rich phenomenology of astrophysical systems including 
active galactic nuclei, X-ray binaries and gamma-ray bursts. Since the radiated en- 
ergy originates from the vicinity of the BH, a fully general-relativistic treatment is 
essential for the modelling of these objects and the flows of plasma in their vicinity. 
Depending on the mass accretion rate, a given system can be found in various 
spectral states, with different radiation mechanisms dominating and varying degrees 
of coupling between radiation and gas |1, 2]. Some supermassive BHs, including 
the primary targets of observations by the Event-Horizon-Telescope Collaboration 
(EHTC")), i.e., Sgr A* and M87, are accreting well below the Eddington accretion 
rate [3, 4]. In this regime, the accretion flow advects most of the viscously released 
energy into the BH rather than radiating it to infinity. Such optically thin, ra- 
diatively inefficient and geometrically thick flows are termed advection-dominated 
accretion flows (ADAFs, see [5-8]) and can be modelled without radiation feedback. 
Next to the ADAF, two additional radiatively inefficient accretion flows (RIAFs) 
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exist: The advection-dominated inflow-outflow solution (ADIOS) [9, 10] and the 
convection-dominated accretion flow (CDAF) [11, 12], which include respectively, 
the physical effects of outflows and convection. Analytical and semi-analytical ap- 
proaches are reasonably successful in reproducing the main features in the spectra 
of ADAFs [see, e.g., 13]. However, numerical general-relativistic magnetohydrody- 
namic (GRMHD) simulations are essential to gain an understanding of the detailed 
physical processes at play in the Galactic Centre and other low-luminosity compact 
objects. 

Modern BH accretion-disk theory suggests that angular momentum transport is 
due to MHD turbulence driven by the magnetorotational instability (MRI) within 
a differentially rotating disk [14, 15]. Recent non-radiative GRMHD simulations of 
BH accretion systems in an ADAF regime have resolved these processes and reveal 
a flow structure that can be decomposed into a disk, a corona, a disk-wind and a 
highly magnetized polar funnel [see, e.g., 16-19]. The simulations show complex 
time-dependent behaviour in the disk, corona and wind. Depending on BH spin, 
the polar regions of the flow contain a nearly force-free, Poynting—flux-dominated 
jet [see, e.g., 17, 18, 20, 21]. 

In addition to having to deal with highly nonlinear dynamics that spans a large 
range in plasma parameters, the numerical simulations also need to follow phenom- 
ena that occur across multiple physical scales. For example, in the MHD paradigm, 
jet acceleration is an intrinsically inefficient process that requires a few thousand 
gravitational radii to reach equipartition of the energy fluxes [22, 23] (purely hy- 
drodynamical mechanisms can however be far more efficient [24]). Jet-environment 
interactions like the prominent HST-1 feature of the radio-galaxy M87 [25-27] oc- 
cur on scales of ~ 5 x 10° gravitational radii. Hence, for a self-consistent picture 
of accretion and ejection, jet formation and recollimation due to interaction with 
the environment [see, e.g., 28], numerical simulations must capture horizon-scale 
processes, as well as parsec-scale interactions with an overall spatial dynamic range 
of ~ 10°. The computational cost of such large-scale grid-based simulations quickly 
becomes prohibitive. Adaptive mesh refinement (AMR) techniques promise an ef- 
fective solution for problems where it is necessary to resolve small and large scale 
dynamics simultaneously. 

Another challenging scenario is presented by radiatively efficient geometrically 
thin accretion disks that mandate extreme resolution in the equatorial plane in order 
to resolve the growth of MRI instabilities. Typically this is dealt with by means of 
stretched grids that concentrate resolution where needed [29, 30]. However, when 
the disk is additionally tilted with respect to the spin axis of the BH [31, 32], lack of 
symmetry forbids such an approach. Here an adaptive mesh that follows the warping 
dynamics of the disk can be of great value. The list of scenarios where AMR can have 
transformative qualities due to the lack of symmetries goes on, the modelling of star- 
disk interactions [33], star-jet interactions [34], tidal disruption events [35], complex 
shock geometries [36, 37], and intermittency in driven-turbulence phenomena [38, 
39], will benefit greatly from adaptive mesh refinement. 

Over the past few years, the development of general-relativistic numerical codes 
employing the 3 + 1 decomposition of spacetime and conservative “Godunov” 
schemes based on approximate Riemann solvers [40—42] have led to great advances 
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in numerical relativity. Many general-relativistic hydrodynamic (HD) and MHD 
codes have been developed [37, 43-59] and applied to study a variety of problems 
in high-energy astrophysics. Some of these implementations provide additional ca- 
pabilities that incorporate approximate radiation transfer [see, e.g., 60-62] and/or 
non-ideal MHD processes [see, e.g., 63, 64]. Although these codes have been applied 
to many astrophysical scenarios involving compact objects and matter [for recent 
reviews see, e.g., 42, 65], full AMR is still not commonly utilised and exploited 
[with the exception of 49, 58, 66]. BHAC attempts to fill this gap by providing a 
fully-adaptive multidimensional GRMHD framework that features state-of-the-art 
numerical schemes. 

Qualitative aspects of BH accretion simulations are code-independent [see, 
e.g., 16, 46, 49], but quantitative variations raise questions regarding numerical 
convergence of the observables [58, 67]. In preparation for the upcoming EHTC ob- 
servations, a large international effort, whose European contribution is represented 
in part by the BlackHoleCam project! [68], is concerned with forward modelling 
of the future event horizon-scale interferometric observations of Sgr A* and M87 at 
submillimeter (EHTC; [69]) and near-infrared wavelengths (VLTI GRAVITY; [70]). 
To this end, GRMHD simulations have been coupled to general-relativistic radiative 
transfer (GRRT) calculations [see, e.g., 71-76]. In order to assess the credibility of 
these radiative models, it is necessary to assess the quantitative convergence of the 
underlying GRMHD simulations. In order to demonstrate the utility of BHAC for the 
EHTC science-case, we therefore validate the results obtained with BHAC against the 
HARM3D code [46, 77] and investigate the convergence of the GRMHD simulations 
and resulting observables obtained with the GRRT post-processing code BHOSS [78]. 

The structure of the paper is as follows. In Sect. 2 we describe the governing 
equations and numerical methods. In Sect. 3 we show numerical tests in special- 
relativistic and general-relativistic MHD. In Sect. 4 the results of 2D and 3D 
GRMHD simulations of magnetised accreting tori are presented. In Sect. 5 we briefly 
describe the GRRT post-processing calculation and the resulting image maps from 
the magnetised torus simulation shown in Sect. 4. In Sect. 6 we present our conclu- 
sions and outlook. 


Throughout this paper, we adopt units where the speed of light, c = 1, the 
gravitational constant, G = 1, and the gas mass is normalised to the central compact 
object mass. Greek indices run over space and time, i.e., (0,1,2,3), and Roman 
indices run over space only, i.e., (1, 2,3). We assume a (—, +, +, +) signature for the 
spacetime metric. Self-gravity arising from the gas is neglected. 


2 GRMHD formulation and numerical methods 
In this section we briefly describe the covariant GRMHD equations, introduce the 
notation used throughout this paper, and present the numerical approach taken in 
our solution of the GRMHD system. The computational infrastructure underlying 
BHAC is the versatile open-source MPI-AMRVAC toolkit [79, 80]. 

In-depth derivations of the covariant fluid- and magneto-fluid dynamical equations 
can be found in the textbooks by [40, 81, 82]. We follow closely the derivation of 


http: //www.blackholecam.org 


Page 3 of 62 


Porth et al. 


the GRMHD equations by [52]. This is very similar to the “Valencia formulation”, 
cf. [40] and [50]. The general considerations of the “3+1” split of spacetime are 
discussed in greater detail in [83-85]. 

We start from the usual set of MHD equations in covariant notation 


V (pu) = 0, 
VaT" =0, (1) 
VF” =0, 


which respectively constitute mass conservation, conservation of the energy- 
momentum tensor 7”, and the homogeneous Faraday’s law. The Faraday tensor 
FY” may be constructed from the electric and magnetic fields E“, B% as measured 
in a generic frame U® as 


F" = UHE” — U” E! — (gh OU Bs, (2) 


where n>? is the fully-antisymmetric symbol (see, e.g., [40]) and g the determinant 
of the spacetime four-metric. The dual Faraday tensor *F¥Y := 4 (—g)7 12n ^ Fy 
is then 


“FAY = UHB” — U” BY — (—g) Vn U, Es . (3) 


We are interested only in the ideal MHD limit of vanishing electric fields in the fluid 
frame u”, hence 


F"”u, =0, (4) 
such that the inhomogeneous Faraday’s law is not required and electric fields are 
dependent functions of velocities and magnetic fields. To eliminate the electric fields 
from the equations it is convenient to introduce vectors in the fluid frame and 
therefore we define the corresponding electric and magnetic field four-vectors as 

e” := Fug, b := “F” u, (5) 
where e” = 0 and we obtain the constraint u,b” = 0. The Faraday tensor is then 


F” = —(—g) 12n > uba, *pHY = btu” — but, (6) 


and we can write the total energy-momentum tensor in terms of the vectors u” and 
b! alone [86] as 


pre Phiotu” u” + Protg”” — bb. (7) 
Here the total pressure Prot = p + 67/2 was introduced, as well as the total specific 


enthalpy htot = h + b?/p. In addition, we define the scalar b? := b”b,, denoting the 
square of the fluid frame magnetic field strength as b? = B? — E?. 
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2.1 3+1 split of spacetime 

We proceed to split spacetime into 3+1 components by introducing a foliation into 
space-like hyper-surfaces 4, defined as iso-surfaces of a scalar time function t. This 
leads to the timelike unit vector normal to the slices 4; [40, 85] 


ny = —aVyt, (8) 


where a is the so-called lapse-function. The four-velocity n” defines the frame of 
the Eulerian observer. If gy, is the metric associated with the four-dimensional 
manifold, we can define the metric associated with each timelike slice as 


Yuv i= Juv + Npn. (9) 
This also allows us to introduce the spatial projection operator 

yE = OF + nl ny (10) 
such that yn, = 0 and through which we can project any four-vector V” (or 
tensor) into its temporal and spatial components. 


Introducing a coordinate system adapted to the foliation X+, the line element is 
given in 3+1 form [87] as 


ds? = —a° dt? + yi;(dax' + B'dt)(da? + pidt), (11) 


where the spatial vector 8" is called the shift vector. Written in terms of coordinates, 
it describes the motion of coordinate lines as seen by an Eulerian observer 


Tipat = zi = p (t, x) dt . (12) 
More explicitly, we write the metric guy and its inverse g’” as 


z —a? + Grp" Bi w af 71/2 Bt Ja? 
Juv B; “ij , g BI Ja? aid _ Bi BI Ja? . 


(13) 


From (13) we find the following useful relation between the determinants of the 
3-metric and 4-metric 


(-g)'? = aq. (14) 


In a coordinate system specified by (11), the four-velocity of the Eulerian observer 
reads 


ny, = (—a,0;), n” = (1/a,—6/a). (15) 
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It is easy to verify that this normalised vector is indeed orthogonal to any space- 
like vector on the foliation X+. Given a fluid element with four-velocity u”, the 
Lorentz factor with respect to the Eulerian observer is?! T := —ulny, = au? and 
we introduce the three-vectors 

i et w a 


US =—+ vi := VV) = i (16) 
T T a’ j p~ 


which denote the fluid three-velocity. 
In the following, purely spatial vectors (e.g., v? = 0) are denoted by Roman 
indices. Note that T = (1 — v?)~1/? with v? = u,v" just as in special relativity. 
Further useful three-vectors are the electric and magnetic fields in the Eulerian 


frame 
E’ := F”n, =a F”, B? := *F”n, = o F” , (17) 


which differ by a factor a from the definitions used in [46, 88]. Writing the general 
Faraday tensor (2) in terms of quantities in the Eulerian frame, the ideal MHD 
condition (4) leads to the well known relation 


EP = Pn" Bou, , (18) 


or put simply: E = B x v (here nij is the standard Levi-Civita antisymmetric 
symbol). Combining (6) with (17), one obtains the transformation between b“ and 
BY as 


pa abut 


bp p? = T(Btv;) 


T i a 


(19) 


which enables the dual Faraday tensor (6) to be expressed in terms of the Eulerian 
fields 
Bru” — Bu” 
fo . 20 
= (20) 
Equation (1) with the Faraday tensor in the form (20) yields the final evolution 
equation for B. The time component of this leads to the constraint ði yB =0 
or put more simply: V - B = 0. Following (19) we obtain the scalar b? as 


3 _ B? + a? (b0)? _ B? 


p T? T? 


} (Bivi)? , (21) 


where B? := B! B;. 

Using the spatial projection operator, the GRMHD equations (1) can be decom- 
posed into spatial and temporal components. We skip ahead over the involved al- 
gebra [see e.g., 52] and directly state the final conservation laws 


O(V7U) + (V7 F") = V7S, (22) 


BlThis quantity is often indicated as w [40, 50] 
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with the conserved variables U and fluxes F" defined as 


D VİD 

p= e i (23) 
T a(S — v D) — B'r 
Bi V' BI — BiVI 


where we define the transport velocity V? := avt — pt. Hence we solve for conservation 
of quantities in the Eulerian frame: the density D := —pu”n,, the covariant three- 
momentum S;, the rescaled energy density 7 = U — D [4] (where U denotes the 
total energy density as seen by the Eulerian observer), and the Eulerian magnetic 
three-fields BI. The conserved energy density U is given by 


U := T” nny = phr? — p+ (E? + B?) (24) 


1 
2 
1 
= phT? — p+ = 
p prs 


[B?(1 +0?) — (Blv;)?] . (25) 


The purely spatial variant of the stress-energy tensor W“) was introduced for ex- 
ample in (23). It reads just as in special relativity 


WY = aT = ph? vu) — EEI — BB? 4 [p+ JE +54) vy) (26) 


B’ BI 


= Stu! poy — a (B B. (27) 


Correspondingly, the covariant three-momentum density in the Eulerian frame is 


Si = Yë n” Tau = phr?vi + Nijr y EÍ BP (28) 
= phT?y; + B?v; — (Bivi) Bi, (29) 


as usual. For the sources S we employ the convenient Valencia formulation without 
Christoffel symbols, yielding 


0 
3aW Ojik + S0; p" —Uðja 
EWY Bi O;vin + W706 — SIOj0 
0 


(30) 


which is valid for stationary spacetimes that are considered for the remainder of this 
work (Cowling approximation). Following the definitions (23) and (30), all vectors 
and tensors are now specified through their purely spatial variants and thus apart 
from the occurrence of the lapse function a and the shift vector 8", the equations 
take on a form identical to the special-relativistic MHD (SRMHD) equations. This 


Using + = U—D instead of u improves accuracy in regions of low energy and enables 
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fact allows for a straightforward transformation from the SRMHD physics module 
of MPI-AMRVAC into a full GRMHD code. 

In addition to the set of conserved variables U, knowledge of the primitive vari- 
ables P(U) is required for the calculation of fluxes and source terms. They are 
given by 


While the transformation U(P) is straightforward, the inversion P(U) is a non- 
trivial matter which will be discussed further in Sect 2.10. Note that just like in 
MPI-AMRVAC , we do not store the primitive variables P but extend the conserved 
variables by the set of auxiliary variables 


A=II,é], (32) 


where € := I'?ph. Knowledge of A allows for quick transformation of P(U). The 
issue of inversion then becomes a matter of finding an A consistent with both P 
and U. 


2.2 Finite volume formulation 
Since BHAC solves the equations in a finite volume formulation, we take the integral 
of equation (22) over the spatial element of each cell f dxtdx?dx* 


facrPuyactastde + fay? detdatde’ = [rPsactactas’ 
(33) 


This can be written (cf. [89]) as 


a,(U AV) + | yl? Fl dada’ — | hl? Fl a 
OV (a14+Aa! /2) OV (a1—Aa1 /2) 
+f yl? Fda da? -f y? F?dz!dr’ 
OV (x? +Ax?/2) OV (x2?— Az? /2) 
+f yt? F8 da! da? -f y? F’dz!dx? 
OV (x3+Ax3/2) OV (a3 — Az? /2) 
= SAV, 
(34) 
with the volume averages defined as 
oT f Uda! dx?dx? 5- J~ P Sdz'dx*dx3 (35) 


AV , AV f 


and 


AV = [rac acter’. (36) 
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We next define also the “surfaces” A.S* and corresponding surface-averaged fluxes 
8 8 


AS$v(ritAzi/2) = O Pdt, (37) 
ƏV (zi+Azi/2) 
and 
=) Sov (wi+a0%/2) per dali 
Fay (ait Aci/2) = Agi ; (38) 


Considering that AV is assumed constant in time, this leads to the evolution equa- 


tion 


rT 1 pl 1 pl il 
OU = -Ay F AS vea T AS lare aar 
p2 2 712 2 (39) 
FAS reaa T AS E aa 
m3 3 m3 3 Q 
FAS laveana T AS lav(es-Az3/2)] +S. 


We aim to achieve second-order accuracy and represent the interface-averaged flux, 
e.g., F'oy(ai4Ar1/2) with the value at the midpoint, change to an intuitive index 
notation Ft; /2j,k) and then arrive at a semi-discrete equation for the average 
state in the cell (i, j, k) as 


dÜ; j,k _ 1 1 1 1 1 
u A FAS aae EAS liiat 
2A 92 2A o2 (40) 
F°AS bami T AS amet 
3A Q3 3A o3 
Fras laksi -FAS laaa + Sij. 


Here the source term S; j,k is also evaluated at the cell barycenter to second-order 
accuracy [90]. Barycenter coordinates 7f are straightforwardly defined as 


i _ [y rdr! ard 
g= AV . 


(41) 


This finite volume form is readily solved with the MPI-AMRVAC toolkit. For ease of 
implementation, we pre-compute all static integrals yielding cell volumes AV, Sur- 
faces AS* and barycenter coordinates. The integrations are performed numerically 
at the phase of initialisation using a fourth-order Simpson’s rule. 

For the temporal update, we interpret the semi-discrete form (40) as an ordinary 
differential equation in time for each cell and employ a multi-step Runge-Kutta 
scheme to evolve the average state in the cell U i,j,k. a procedure also known as 
“method of lines”. At each sub-step, the point-wise interface fluxes F* are obtained 
by performing a limited reconstruction operation of the cell-averaged state U to the 
interfaces (see Sect. 2.8) and employing approximate Riemann solvers, e.g., HLL or 
TVDLF (Sect. 2.9). 
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Several temporal update schemes are available: simple predictor-corrector, third- 
order Runge-Kutta (RK) RK3 [91] and the strong-stability preserving s-step, pth- 
order RK schemes SSPRK(s,p) schemes: SSPRK(4,3), SSPRK(5,4) due to [92].5] 


2.3 Metric data-structure 

The metric data-structure is built to be optimal in terms of storage while remaining 
convenient to use. Since the metric and its derivatives are often sparsely populated, 
the data is ultimately stored using index lists. For example, each element in the 
index list for the four-metric gy holds the indices of the non-zero element together 
with a Fortran90 array of the corresponding metric coefficient for the grid block. A 
summation over indices, e.g., “lowering” can then be cast as a loop over entries in the 
index-list only. For convenience, all elements can also be accessed directly over intu- 
itive identifiers which point to the storage in the index list, e.g., mfg (mu, nu) Zelem 
yields the grid array of the g,,, metric coefficients as expected. Similarly, the lower- 
triangular indices point to the transposed indices in the presence of symmetries. In 
addition, one block of zeros is allocated in the metric data-structure and all zero 
elements are set to point towards it. An overview of the available identifiers is given 
in Table 1. 


Table 1 Elements of the metric data-structure 


Symbol | Identifier Index list 

Juv mg (mu ,nu) m/nnonzero, mnonzero(inonzero) 

a m%alpha - 

Bt m%beta(i) m%⁄nnonzeroBeta, m/nonzeroBeta(inonzero) 

VY m/sqrtgamma - 

yd mAgammainv(i,j) | - 

Bi m%betaD (i) - 

On Vig m/dgdk (i,j,k) m/nnonzeroDgDk, mAnonzeroDgDk (inonzero) 

0; B m/DbetaiDj(i,j) | ménnonzeroDbetaiDj, m/nonzeroDbetaiDj (inonzero) 
ðja m4DalphaDj (j) m/nnonzeroDalphaDj, mnonzeroDalphaDj (inonzero) 
0 mýzero - 


As a consequence, only 14 grid functions are required for the Schwarzschild co- 
ordinates and 29 grid functions need to be allocated in the Kerr-Schild (KS) case. 
This is still less than half of the 68 grid functions which a brute-force approach 
would yield. The need for efficient storage management becomes apparent when we 
consider that the metric is required in the barycenter as well as on the interfaces, 
thus multiplying the required grid functions by a factor of four for three-dimensional 
simulations (yielding 116 grid functions in the KS case). 

In order to eliminate the error-prone process of implementing complicated func- 
tions for metric derivatives, BHAC can obtain derivatives by means of an accurate 
complex-step numerical differentiation [93]. This elegant method takes advantage 
of the Cauchy-Riemann differential equations for complex derivatives and achieves 
full double-precision accuracy, thereby avoiding the stepsize dilemma of common 
finite-differencing formulae [94]. The small price to pay is that at the initialisation 
stage, metric elements are provided via functions of the complexified coordinates. 
However, the intrinsic complex arithmetic of Fortran90 allows for seamless imple- 


mentation. 


lFor implementation details, see [80]. 
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To promote full flexibility in the spacetime, we always calculate the inverse metric 
7J using the standard LU decomposition technique [95]. As a result, GRMHD sim- 
ulations on any metric can be performed after providing only the non-zero elements 
of the three-metric yi; (x1, x?, x), the lapse function a(xt, 2”, x°) and the shift vec- 
tor 8*(x2',x?,x3). As an additional convenience, BHAC can calculate the required 


elements and their derivatives entirely from the four-metric guu (x°, £1, x7, x). 


2.4 Equations of state 
For closure of the system (1)-(4), an equation of state (EOS) connecting the specific 
enthalpy h with the remaining thermodynamic variables h(p, p) is required [40]. The 


currently implemented closures are 


e Ideal gas: h(p,p) =1+ = Y P with adiabatic index ĝ. 
g p a 
K(07! 
e Synge gas: h(Q) = a where the relativistic temperature is given by 


© = p/p and K,, denotes the modified Bessel function of the second kind. 
In fact, we use an approximation to the previous expression that does not 
contain Bessel functions [see 79, 96]. 

e Isentropic flow: Assumes an ideal gas with the additional constraint p = Kp’, 
where the pseudo-entropy k may be chosen arbitrarily. This allows one to 
omit the energy equation entirely and only the reduced set P = {p, v’, BÍ} is 
solved. 


As long as h(p,p) is analytic, its implementation in BHAC is straightforward. 


2.5 Divergence cleaning and augmented Faraday's law 

To control the V-B = 0 constraint on AMR grids, we have adopted a con- 
straint dampening approach customarily used in Newtonian MHD [97]. In this ap- 
proach, which is usually referred as Generalized Lagrangian Multiplier (GLM) of 
the Maxwell equations (but is also known as the “divergence-cleaning” approach), 
we extend the usual Faraday tensor by the scalar ¢, such that the homogeneous 
Maxwell equation reads 


VU CEM — bgt”) = —Knko, (42) 


and the scalar ¢ follows from contraction ¢ = (*F#” — dg"”)nynv. Naturally, for 
ġ — 0, the usual set of Maxwell equations is recovered. It is straightforward to show 
[see, e.g., 98] that (42) leads to a telegraph equation for the constraint violation 
parameter ¢ which becomes advected at the speed of light and decays on a timescale 
1/K. With the modification (42), the time-component of Maxwell’s equation now 
becomes an evolution equation for ¢. After some algebra (see Appendix A), we 
obtain 


3T + OL VW(OB" — $8')| =- Vand — 708" 


er (43) 
— zv IB On Vig + YB Oia. 
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Equivalently, the modified evolution equations for Bê (see Appendix B) read 


On (./7B?) + 0; (\/7(V' BY — Vİ B’ — B'B!)) = —\/7B'0;8) — yya did. 
(44) 


Now equation (44) replaces the usual Faraday’s law and (43) is evolved alongside 
the modified MHD system. Due to the term 0;¢ on the right hand side of equation 
(44), the new equation is non-hyperbolic. Hence, numerical stability can be a more 
involved issue than for hyperbolic equations. We find that the numerical stability 


of the system is enhanced when using an upwinded discretisation for 0;¢. Note that 
al 
VY 


Equations (43) and (44) are in agreement with [63] when accounting for 
37’ 0i71m and taking the ideal MHD limit. 


2.6 Flux-interpolated Constrained Transport 

As an alternative to the GLM approach, the V - B = 0 constraint can be enforced 
using a cell-centred version of Flux-interpolated Constrained Transport (FCT) con- 
sistent with the finite volume scheme used to evolve the hydrodynamic variables. 
Constrained Transport (CT) schemes aim to keep to zero at machine precision the 
sum of the magnetic fluxes through all surfaces bounding a cell, and therefore (in 
the continuous limit) the divergence of the magnetic field inside the cell. In the 
original version [99] this is achieved by evolving the magnetic flux through the cell 
faces and computing the circulation of the electric field along the edges bounding 
each face. Since each edge appears with opposite signs in the time update of two 
faces belonging to the same cell, the total magnetic flux leaving each cell is con- 
served during evolution. The magnetic field components at cell centers, necessary 
for performing the transformation from primitive to conserved variables and vice- 
versa, are then found using interpolation from the cell faces. [100] showed that it is 
possible to find cell centred variants of CT schemes that go from the average field 
components at the cell center at a given time to those one (partial) time step ahead 
in a single step, without the need to compute magnetic fluxes at cell faces. The 
CT variant known as FCT is particularly well suited for finite volume conservative 
schemes as that employed by BHAC, as it calculates the electric fields necessary for 
the update as an average of the fluxes given by the Riemann solver. In this way, 
the time update for its cell centred version can be written using a form similar to 
(40). For example, for the update of the Bt component, we obtain 


dB}. 1 
i,j,k * 2 2 
= F*-AS 
dt AV; j,k | 


* 2 2 
i,g+1/2,k — F* "AS eee? 


*3 3 * 3 3 
F*”AS hoe E AS ere y 


where the modified fluxes in the x!-direction are zero and the remaining fluxes are 


calculated as 
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The derivation of equations (45) and (46) from the staggered version with magnetic 
fields located at cell faces is given in Appendix C. Since magnetic fields are stored 
at the cell center and not at the faces, the divergence conserved by the FCT method 
corresponds to a particular discretisation 


savy - B) 


i+1/2,j+1/2,k+1/2 


B!A B?A B3A 
`». [ ij i oe {jy R A is = a 
l ,l2,l3=0,1 s < T Jitlj+la,k+ls 
(47) 
where 
AV" |i41/2,j+1/2,k+1/2 = 5 AV |i 41; jtlo,k-tls : (48) 


11 ,l2,13=0,1 


Equation (47) is closely related to the integral over the surface of a volume contain- 
ing eight cells in 3D (see Appendix D for the derivation), and it reduces to equation 
(27) from [100] in the special case of Cartesian coordinates. As mentioned before, 
this scheme can maintain V - B = 0 to machine precision only if it was already 
zero at the initial condition. The corresponding curl operator used to setup initial 
conditions is derived in Appendix D. 

In its current form, BHAC cannot handle both constrained transport and AMR. 
The reason is that special prolongation and restriction operators are required in or- 
der to avoid the creation of divergence when refining or coarsening. Due to the lack 
of information about the magnetic flux on cell faces, the problem of finding such 
divergence-preserving prolongation operators becomes underdetermined. However, 
storing the face-allocated (staggered) magnetic fluxes and applying the appropriate 
prolongation and restriction operators requires a large change in the code infras- 
tructure on which we will report in an accompanying work. 


2.7 Coordinates 

Since one of the main motivations for the development of the BHAC code is to 
simulate BH accretion in arbitrary metric theories of gravity, the coordinates and 
metric data-structures have been designed to allow for maximum flexibility and 
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can easily be extended. A list of the currently available coordinate systems is given 
in Table 2. In addition to the identifiers used in the code, the table lists whether 
numerical derivatives are used and whether the coordinates are initialised from the 
three-metric or the four-metric. The less well-known spacetimes and coordinates are 


described in the following subsection. 


Table 2 Coordinates available in BHAC 


Coordinates Identifier | Num. derivatives Init. guv 
Cartesian cart No No 
Boyer-Lindquist bl No No 
Kerr-Schild ks No No 
Modified Kerr-Schild mks No No 
Cartesian Kerr-Schild cks Yes Yes 
Rezzolla & Zhidenko parametrization [101] rz Yes No 
Horizon penetrating Rezzolla & Zhidenko coordinates | rzks Yes Yes 
Hartle-Thorne [102] ht Yes Yes 


2.7.1 Modified Kerr-Schild coordinates 
Modified KS coordinates were introduced by e.g., [17] with the purpose of stretching 
the grid radially and being able to concentrate resolution in the equatorial region. 


The original coordinate transformation is equivalent to: 


rxs(s) = Ro + e°, (49) 


Oxs(V) = V + sin(2v) , (50) 


where Ro and h are parameters which control, respectively, how much resolution is 
concentrated near the horizon and near the equator. 

Unfortunately, the inverse of J(@) is a transcendental equation that has to be 
solved numerically. To avoid this complication and still capture the functionality of 
the modified coordinates, we instead use the following 0— transformation 


2hd 
Oxs(0) =V + 


(m — 20)(r — 0). (51) 


Now the solution to the cubic equation can be expressed in closed-form, and the 
only real root reads 


0(O0xs) = 


1 a3 (  272(3n)?/8(h—1) _ 2?/94/3R(Oxs) 
12" R(Oxs) h 


| ov) ; (52) 


where 


R(Oxs) = [ay=sr [-108402- + 108rhôkgs + (h — 4) (2rh + T)?] 
1/3 
O(n — 2s) . (53) 


This is compared with the original version (50) in Fig. 1 and shows a good match 
between the two versions of modified Kerr-Schild coordinates. The radial back- 
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transformation follows trivially as 
s(rgs) = ln(rgs — Ro), (54) 
and the derivatives for the diagonal Jacobian are 


OsTKS =g" (55) 
sks = 1+2h +12h((8/r)? — 0/n). (56) 


— Orig 


— New 


0.0 i 0.5 1.0 i 1.5 2.0 2.5 3.0 


Figure 1 6-grid stretching functions comparing the transcendental function 0(@xs) (solid red 
curves) with the cubic approach (solid blue curves) for h = 0.9. We also give the respective 
derivatives d0/dv (dashed). 


With these transformations, we obtain the new metric gmks = J'gxgJ. Note 
that whenever the parameters h = 0 and Ro = 0 are set, our MKS coordinates 
reduce to the standard logarithmic Kerr-Schild coordinates. 


2.7.2 Rezzolla & Zhidenko parametrization 

The Rezzolla-Zhidenko parameterisation [101] has been proposed to describe 
spherically-symmetric BH geometries in metric theories of gravity. In particular, 
using a continued-fraction expansion (Padé expansion) along the radial coordinate, 
deviations from general relativity can be expressed using a small number of coeffi- 
cients. The line element reads 


2 
ds? = —N(r)at? + 2 ap? 4 p2a6? + 7? sin? A (57) 
N?(r) 


with N(r) and B(r) being functions of the radial coordinate r. The radial position of 
the event horizon is fixed at r = ro > 0 which implies that N (ro) = 0. Furthermore, 
the radial coordinate is compactified by means of the dimensionless coordinate 


g:=1-—, (58) 


r 


in which x = 0 corresponds to the position of the event horizon, while x = 1 
corresponds to spatial infinity. Through this dimensionless coordinate, the function 
N can be written as 


N? = z A(z), (59) 


Page 15 of 62 


Porth et al. 


where A(x) >0 for 0< x< 1. Introducing additional coefficients €, an, and bn, 
the metric functions A and B are then expressed as follows 


A(x) = 1—e1—2) + (ao —6)(1— 2)? + A(x)(1— 2), (60) 
1+b9(1—2) + B(x)(1—2)?. (61) 


3 
È 
l 


Here A and B are functions describing the metric near the event horizon and at 
spatial infinity. In particular, A and B have rapid convergence properties, that is 
by Padé approximants 


7 ay = by 


A(x) = ~ mnr S) B(x) = = p (62) 
1+ — aT 1+—— 
ea 14 23% 
ben ie ee 
where @1, @2,43,... and by, bo, b3,... are dimensionless coefficients that can, in prin- 


ciple, be constrained from observations. The dimensionless parameter e€ is fixed by 
the ADM mass M and the coordinate of the horizon rp. It measures the deviation 
from the Schwarzschild case as 


2M — 2M 
= Mat _ (: J (63) 


TO TO 


It is easy to see that at spatial infinity (x = 1), all coefficients contribute to (62), 
while at event horizon only the first two terms remain, t.e. 


A(0)=a, B(0)=b. (64) 


Given a number of coefficients, any spherical spacetime can hence directly be simu- 
lated in BHAC. For example, the coefficients in the Rezzolla-Zhidenko parametriza- 
tion for the Johannsen-Psaltis [103] metric and for Einstein-Dilaton BHs [104] have 
already been provided in [101]. Typically, expansion up to ag, be yields sufficient nu- 
merical accuracy for the GRMHD simulations. The first simulations in the related 
horizon penetrating form of the Rezzolla-Zhidenko parametrization are discussed in 
[105]. 


2.8 Available reconstruction schemes 

The second-order finite volume algorithm (40) requires numerical fluxes centered 
on the interface mid-point. As in any Godunov-type scheme [see e.g., 88, 106], 
the fluxes are in fact computed by solving (approximate) Riemann problems at 
the interfaces (see Sect. 2.9). Hence for higher than first-order accuracy, the fluid 
variables need to be reconstructed at the interface by means of an appropriate spa- 
tial interpolation. Our reconstruction strategy is as follows. 1) Compute primitive 
variables P from the averages of the conserved variables U located at the cell 
barycenter. 2) Use the reconstruction formulae to obtain two representations for 
the state at the interface, one with a left-biased reconstruction stencil P” and the 
other with a right-biased stencil P®. 3) Convert the now point-wise values back to 
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their conserved states U” and UF. The latter two states then serve as input for the 
approximate Riemann solver. 

A large variety of reconstruction schemes are available, which can be grouped into 
standard second-order total variation diminishing (TVD) schemes like “minmod”, 
“vanLeer”, “monotonized-central” , “woodward” and “koren” [see 79, for details] and 
higher order methods like the third-order methods “PPM” [107], “LIMO3” [108] 
and the fifth-order monotonicity preserving scheme “MP5” due to [109]. While 
the overall order of the scheme will remain second-order, the higher accuracy of 
the spatial discretisation usually reduces the diffusion of the scheme and improves 
accuracy of the solution [see, e.g., 80]. For typical GRMHD simulations with near- 
evacuated funnel/atmosphere regions, we find the PPM reconstruction scheme to 
be a good compromise between high accuracy and robustness. For simple flows, 
e.g., the stationary toroidal field torus discussed in Sect. 3.4, the compact stencil 
LIMO3 method is recommended. 


2.9 Characteristic speed and approximate Riemann solvers 

The time-update of BHAC proceeds in a dimensionally unsplit manner, thus at each 
Runge-Kutta substep the interface-fluxes in all directions are computed based on the 
previous substep. The state is then advanced to the next substep with the combined 
fluxes of the cell. To compute these fluxes from the reconstructed conserved variables 
at the interface U” and UË, we provide two approximate Riemann solvers: 1 ) the 
Rusanov flux, also known as Total variation diminishing Lax-Friedrichs scheme 
(TVDLF) which is based on the largest absolute value of the characteristic waves 
normal to the interface c’, and 2) the HLL solver [110], which is based on the 
leftmost (c’_) and rightmost (c',) waves of the characteristic fan with respect to 
the interface. The HLL upwind flux function for the conserved variable u € U is 
calculated as 


ra) > So 
F?(u) = FU") ; É <0 (65) 
Fi(U",U"®) ; otherwise 


where 


É, F’ (u") -£ F’ (U®) + a ch (uÈ — u!) 


F(U}, UP) := (66) 
and we set in accordance with [111]: c+ = min (AŁ_,AR_), ci, = max (AG Ae) 

The TVDLF flux is simply 

i Ll pie icp TR lis R L 
Fi(u) = 5 [F (U) + F(U®)| - Se (uP ul) . (67) 

with c’ = max (| |, |é). 

In addition to these two standard approximate Riemann solvers, we also provide 
a modified TVDLF solver that preserves positivity of the conserved density D. The 
algorithm was first described in the context of Newtonian hydrodynamics by [112] 
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and was successfully applied in GRHD simulations by [113]. It takes advantage of the 
fact that the first-order Lax-Friedrichs flux F*™°(u) is positivity preserving under 
a CFL condition CFL < 1/2. Hence the fluxes can be constructed by combining the 
high order flux F*#°(u) (obtained e.g., by PPM reconstruction) and F®ŁO (u) such 
that the updated density does not fall below a certain threshold.) Specifically, the 
modified fluxes read 


Fi(u) = 6F*H9(u) + (1 — 0) F19 (u), (68) 


where @ € [0,1] is chosen as a maximum value which ensures positivity of the 
cells adjacent to the interface (see [112] for details of its construction). Note that 
although we only stipulate the density be positive, the formula (68) must be applied 
to all conserved variables u € U. 


In relativistic MHD, the exact form of the characteristic wave speeds A+ involves 


solution of a quartic equation [see, e.g., 86] which can add to the computational 
overhead. For simplicity, instead of calculating the exact characteristic velocities, 
we follow the strategy of [46] who propose a simplified dispersion relation for the 
fast MHD wave w? = a?k?. As a trade-off, the simplification can overestimate the 
wavespeed in the fluid frame by up to a factor of 2, yielding a slightly more diffusive 
behaviour. The upper bound a for the fast wavespeed is given by 


2. ae 2 242 
a” = C3 + Ca — CsCa > (69) 


which depends on the usual sound speed and Alfvén speed 


2 b? 


2 aP C, = ————z 
sS Vh’ a ph +b?’ 


c (70) 
here given for an ideal EOS with adiabatic index Ẹ. As pointed out by [52], the 3+1 


structure of the fluxes leads to characteristic waves of the form 


Mo aii- G, (71) 


where į is the characteristic velocity in the corresponding special relativistic sys- 


tem (a > 1, 6f + 0). 
For the simplified isotropic dispersion relation, the characteristics can then be 
obtained just like in special relativistic hydrodynamics [see, e.g., 89, 114, 115] 


(1 — a°) vit y? (1 — v?) la — va?) gë” — (1 — a?) (vi)? 

i : 72 

(1 — v2 a?) ee) 
(In the general-relativistic hydrodynamic WhiskyTHC code [54, 55], this desirable 
property allows to set floors on density close to the limit of floating point precision 


> 
l 


~ 107 ti pref. 
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2.10 Primitive variable recovery 

It is well-known that the nonlinear inversion P(U) is the Achilles heel of any rel- 
ativistic (M)HD code and sophisticated schemes with multiple backup strategies 
have been developed over the years as a consequence (e.g., [116], [117], [77], [118], 
[119], [120]). Here we briefly describe the methods used throughout this work and 
refer to the previously mentioned references for a more detailed discussion. 


2.10.1 Primary inversions 

Two primary inversion strategies are available in BHAC. The first strategy, which we 
denote by “1D”, is a straightforward generalisation of the one-dimensional strategy 
described in [121]. It involves a non-linear root finding algorithm which is imple- 
mented by means of a Newton-Raphson scheme on the auxiliary variable €. Once £ 
is found, the velocity follows from (29) 


Si Bi(BIS;) 


~ (E+B?) (E+ B?)’ 


a 


(73) 


and we calculate the second auxiliary variable T = (1 — v?)~!/? so that p = D/T. 
The thermal pressure p then follows from the particular EOS in use (Sect. 2.4). For 
example, for an ideal EOS we have 


ny a 


For details of the consistency checks and bracketing, we refer the interested reader 
to [121]. 

In addition to the 1D scheme, we have implemented the “2DW” method of [52, 
116]. The 2DW inversion simultaneously solves the non-linear equations (25) and the 
square of the three-momentum $7, following (29) by means of a Newton-Raphson 
scheme on the two variables € and v?. Among all inversions tested by [116], the 
2DW method was reported as the one with the smallest failure rate. We find the 
same trend, but also find that the lead of 2DW over 1D is rather minor in our tests. 

With two distinct inversions that might fail under different circumstances, one can 
act as a backup strategy for the other. Typically we first attempt a 2DW inversion 
and switch to the 1D method when no convergence is found. The next layer of 
backup can be provided by the entropy method as described in the next section. 


2.10.2 Entropy switch 

To deal with highly magnetised regions, [60, 77] introduced the advection of entropy 
to provide a backup strategy for the primitive variable recovery. Similar to [60, 77], 
alongside the usual fluid equations, BHAC can be configured to solve an advection 
equation for the entropy S 


Vp Su” =0, (75) 
where we define 


Spa, (76) 
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given the adiabatic index 7. This leads to the evolution equation 
SWS + O/Y(av' — BTS =0, (77) 


for the conserved quantity TS. The primitive counterpart is the actual entropy 
k = p/p’, which can be recovered via x = 'S/D. In case of failure of the primary 
inversion scheme, using the advected entropy «, we can attempt a recovery of prim- 
itive variables which does not depend on the conserved energy. Note that after the 
primitive variables are recovered from the entropy, we need to discard the conserved 
energy and set it to the value consistent with the entropy. On the other hand, after 
each successful recovery of primitive variables, the entropy is updated to K = p/p’, 
which is then advected to the next step. In addition, entropy-based inversion can 
be activated whenever 8 = 2p/b? < 107? since the primary inversion scheme is 
likely to fail in these highly magnetised regions. Tests of the dynamic switching of 
the evolutionary equations are described in Sect. 3.3. In GRMHD simulations of 
BH accretion, the “entropy region” is typically located in the BH magnetosphere, 
which is strongly magnetised and the error due to missing shock dissipation is thus 
expected to be small. 

In the rare instances where the entropy inversion also fails to converge to a physical 
solution, the code is normally stopped. To force a continuation of the simulation, 
last resort measures that depend on the physical scenario can be employed. Often 
the simulation can be continued when the faulty cell is replaced with averages of 
the primitive variables of the neighbouring healthy cells as described in [79]. In the 
GRMHD accretion simulations described below, failures could happen occasionally 
in the highly magnetised evacuated “funnel” region close to the outer horizon where 
the floors are frequently applied. We found that the best strategy is then to replace 
the faulty density and pressure values with the floor values and set the Eulerian 
velocity to zero. Note that in order to avoid generating spurious V - B, the last 
resort measures should never modify the magnetic fields of the simulation. 


2.11 Adaptive Mesh refinement 

The computational grid employed in BHAC is provided by the MPI-AMRVAC toolkit 
and constitutes a fully adaptive block based (oct-) tree with a fixed refinement factor 
of two between successive levels. That is, the domain is first split into a number of 
blocks with equal amount of cells (e.g., 103 computational cells per block). Each 
block can be refined into two (1D), four (2D) or eight (3D) child-blocks with an 
again fixed number of cells. This process of refinement can be repeated ad libitum 
and the data-structure can be thought of a forest (collection of trees). All operations 
on the grid, for example time-update, IO and problem initialisation are scheduled 
via a loop over a space-filling curve. We adopt the Morton Z-order curve for ease of 
implementation via a simple recursive algorithm. 

Currently, all cells are updated with the same global time-step and hence load- 
balancing is achieved by cutting the space-filling curve into equal sections that 
are then distributed over the MPI-processes. The AMR strategy just described is 
applied in various astrophysical codes, for example codes employing the PARAMESH li- 
brary [122-124], or the recent Athena++ framework [see, e.g., 125]. Compared to a 
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patch-based approach [see, e.g., 126], the block based AMR has several advantages: 
1) well-defined boundaries between neighbouring grids on different levels ,2) data 
is uniquely stored and updated, thus no unnecessary interpolations are performed, 
and 3) simple data-structure, e.g., straightforward integer arithmetic can be used to 
locate a particular computational block. For in-depth implementation details such 
as refinement /prolongation operations, indexing and ghost-cell exchange, we refer to 
[79]. Prolongation and restriction can be used on conservative variables or primitive 
variables. Typically primitive variables are chosen to avoid unphysical states which 
can otherwise result from the interpolations in conserved variables. The refinement 
criteria usually adapted is the Léhner’s error estimator [127] on physical variables. 
It is a modified second derivative, normalised by the average of the gradient over 
one computational cell. The multidimensional generalization is given by 
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(78) 


The indices p,q run over all dimensions p,q = 1,...,Np. The last term in the 
denominator acts as a filter to prevent refinement of small ripples, where fwave is 
typically chosen of order 107?. This method is also used in other AMR codes such 
as FLASH [128], RAM [124], PLUTO [126] and ECHO [66]. 


3 Numerical tests 
3.1 Shock tube test with gauge effect 
The first code test is considered in flat spacetime and therefore no metric source 
terms are involved. Herein we perform one-dimensional MHD shock tube tests with 
gauge effects by considering gauge transformations of the spacetime. Shock tube 
tests are well-known tests for code validation and emphasise the nonlinear behaviour 
of the equations, as well as the ability to resolve discontinuities in the solutions [see, 
e.g., 50, 52]. 

The initial condition is given as 


(1,1, 0.5, 1) x<0, 


79 
(0.125, 0.1,0.5,-1) «>0, 2) 


(p, p, B”, B”) “| 


and all other quantities are zero. In order to check whether the covariant fluxes are 
correctly implemented, we use different settings for the flat spacetime as detailed 
in Table 3. 


Table 3 Shock tube with gaugeeffect setups. 


Case a BY Yua V22 V33 
A i (000) 1 í 1 
B 2 (000) 1 #21 1 
C 1 (0400) 1 #1 1 
D 1 (000) 4 #1 1 
E 1 (000) 14 4 1 
F 2 (0400) 4 9 1 
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In the simulations, an ideal gas EOS is employed with an adiabatic index of 7 = 2. 
The 1D problem is run on a uniform grid in «—direction using 1024 cells spanning 
over x € [—1/2,1/2]. The simulations are terminated at t = 0.4. For the spatial 
reconstruction, we adopt the second order TVD limiter due to Koren [129]. Fur- 
thermore, RK3 timeintegration is used with Courant number set to 0.4. 

Case A is the reference solution without modification of fluxes due to the three- 
metric, lapse or shift!”], By means of simple transformations of flat-spacetime, all 
other cases can be matched with the reference solution. Case B will coincide with 
solution A if B is viewed at t/2 = 0.2. Case C will agree with case A when it is 
shifted in positive x—direction by da = 6*t = 0.16. For case D, we rescale the 
domain as x € [—1/4,1/4] and initialise the contravariant vectors as B’? = B® /2. 
The state at t = 0.4 should agree with case A when the domain is multiplied by the 
scale factor hy = 2. For case E we initialise B’” = BY/2 and case F is initialised 
similarly as B’? = B*/2, BY = BY/3. 

In general, all cases agree very well with the rescaled solution. To give an example, 
Fig. 2 shows the rescaled simulation results of case F compared to the reference 
solution of case A. This test demonstrates the shock-capturing ability of the MHD 
code and enables us to conclude that the calculation of the covariant fluxes has 


been implemented correctly. 


1L 1 1 f fi -0.8L f f f f f f f L n 
—0.4 -0.2 0.0 02 0.4 —0.4 -0.2 0.0 02 04 -0.4 -0.2 00 02 04 


Figure 2 1D plots of density p, gas pressure p, Lorentz factor I’, velocity components v” and vY, 
and the y-component of the magnetic field for the shock tube test at t = 0.4. The reference 
solution of case A is shown as a solid black line and the rescaled solution of case F is overplotted 
as red squares. 


3.2 Boosted loop advection 

In order to test the implementation of the GLM-GRMHD system, we perform the 
advection of a force-free flux-tube with poloidal and toroidal components of the 
magnetic field in a flat spacetime. 


We note that for the reference solution we have relied here on the extensive set 
of tests performed in flat spacetime within the mp1-amavac framework; however, we 
could also have employed as reference solution the “exact” solution as derived in 
Ref. [130]. 
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The initial equilibrium configuration of a force-free flux-tube is given by a mod- 
ified Lundquist tube [see e.g., 131], where we avoid sign changes of the vertical 
field component B* with the additive constant C = 0.01. Pressure and density are 
initialized as constant throughout the simulation domain. The initial pressure value 
is obtained from the central plasma-beta Bo = B?(0)/2p, where B is the magnetic 
field in the co-moving system. The density is set to p = p/2 yielding a relativistic 
hot plasma. Consequently, an adiabatic index y = 4/3 is used. We set Bọ = 0.01, 
which results in a high magnetisation co = B?(0)/(pc? + 4p) ~ 25. The equations 
for the magnetic field for r < 1 read 


Be(r) = Ji (œr), (80) 


B*(r) = y Jalar +C, (81) 
and 


B?(r)=0, (82) 
B? (r) = y Jola)? +C, (83) 


otherwise, where Jo and J; are Bessel functions of zeroth and first order respectively 
and the constant a; ~ 3.8317 is the first root of Jo. 

This configuration is then boosted to the frame moving at velocity v = 
V2(—ve, —ve, 0) and we test values of ve between 0.5¢ and 0.99c. 


Standard Lorentz transformation rules result in 


T2 
r+1 


r=r' + (T-1)(r'.-njn-Ttvn, B'=TB-— B(B-B), (84) 
where t’ can be set to zero and where we assumed a vanishing electric field in 
the co-moving system. Therefore relativistic length contraction gives the loop a 
squeezed appearance. A simulation domain (x,y) € [—1,1] at a base resolution of 
Nz x Ny = 64? is initialised with an additional three levels of refinement. We advect 
the loop for one period (P = 2\/2/v) across the domain where periodic boundary 
conditions are used. 

The advection over the coordinates can be counteracted by setting the shift vector 
appropriately, i.e. 8 = —v. This is an important consistency check of the imple- 
mentation. Figure 3 shows the initial and final states of the force-free magnetic 
flux-tube advected for one period and for the case with spacetime shifted against 
the advection velocity. The advected and counter-shifted cases are in good agree- 
ment, with only the truly advected case being slightly more diffused, the effect of 
which is reflected in the activation of more blocks on the third AMR level. 

To investigate the numerical accuracy the Lı and ZL. norms of the out-of-plane 
magnetic field component B7, as well as the divergence of magnetic field between the 
initial state and the simulation at a time after one advection period with different 


resolutions as seen in Fig.4 are checked. The error norms from analytically known 
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Be, t=0 Be , t=5.65685 


1.0 1.0 
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B? , t=5.65685 


Figure 3 Force-free magnetic loop with diagonal boost velocity |v| = 0.5c. Top: No shift, the loop 
is advected for one period. Bottom: The shift vector just opposes the (diagonal) advection 
velocity, |v| = 0.5, hence the loop remains stationary with respect to the grid. Base resolution is 
64? cells with a total of three grid levels. The color shows the strength of the out-of-plane field 
component B? and white lines are in-plane field lines of (B®, B7). Blocks containing 8? cells are 
indicated. 
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solutions u* are defined as 


1 1 
AE tiep — ——— u“ ,/ydx'dz’dz?| , 85 
1( Neells 2 ihk AVi j,k Vi,j,k v7 l l 
1 _ 
Doo (ts) = maxis e tink — Ry [ u* /yda'da*da*) , (86) 
Ud; i,j,k 


where the summation, respectively maximum operation, includes all cells in the 
domain and the integrals are performed over the volume of the cell AV; j,- In 
this sense, the reported errors correspond to the mean and maximal error in the 
computational domain. We should note that for the test of convergence, we use 
a uniform grid and choose v = 0.5V0.5(1,1,0),8 = V0.5(1,0,0) resulting in an 
advection in direction of the upper-left diagonal. A TVD “Koren” limiter is chosen. 
As expected, the convergence is second order for all cases. 


Figure 4 Error of the out-of-plane magnetic field component B* (left) and divergence of B 
(right). For this test, we chose v = 0.5V/0.5(1, 1,0) and 6 = V0.5(1,0,0), resulting in an 
advection in the direction of the upper-left diagonal. 


3.3 Magnetised spherical accretion 
A useful stress test for the conservative algorithm in a general-relativistic setting 
is spherical accretion onto a Schwarzschild BH with a strong radial magnetic field 
[46]. The steady-state solution is known as the Michel accretion solution [132] and 
represents the extension to general relativity of the corresponding Newtonian solu- 
tion by [133]. The steady-state spherical accretion solution in general relativity is 
described in a number of works [see, e.g., 40, 43]. It is easy to show that the solu- 
tion is not affected when a radial magnetic field of the form B” œ y~!/? is added 
[45]. This test challenges the robustness of the code and of the inversion procedure 
P(U) in particular. The calculation of the initial condition follows that outlined in 
[43]. Here, we parametrize the field strength via ø = b?/p at the inner edge of the 
domain (r = 1.9M). The simulation is setup in the equatorial plane using MKS 
coordinates corresponding to a domain of rks € [1.9, 20] M. The analytic solution 
remains fixed at the inner and outer radial boundaries. We run two cases, case 1 
with magnetisation up to o = 10° and case 2 with a very high magnetisation reach- 
ing up to ø = 104. Since the problem is only 1D, the constraint V-B = 0 has a 
unique solution which gets preserved via the FCT algorithm. 

Figure 5 illustrates the profiles for ø = 10° and two inversion strategies: 2DW 
(black +) and 2DW with entropy switching in regions of high magnetization b?/2p > 
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Figure 5 Magnetized Bondi flow at t = 100 M in MKS coordinates with o = 10° at the inner 
edge of the domain. The black solid curve indicates the initial exact solution. We show two 
realisations with resolution N, = 100. Black crosses are with the standard treatment for the 
inversion. Red crosses switch to the entropy evolution at 6 < 10~? (here in the middle of the 
domain). In particular, the error in the radial three-velocity v” decreases when switching to the 
entropy evolution. 
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100 (red x). With the exception of the radial three-velocity near the BH horizon 
(r <5 M), in both cases the simulations maintain well the steady-state solution." 
Comparing theses results with and without entropy switching, the entropy strategy 
actually keeps the solution closer to the steady-state solution (solid black curves) 
even though the change of inversion strategy occurs in the middle of the domain, 


r~ 10. 


Error 


Figure 6 Error of density p in the highly magnetised Bondi flow: ø = 10° (left) and o = 104 
(right). The black data points are obtained with the standard 2D inversion and the red datapoints 
switch to the entropy inversion at 8 < 10~?. One can see that both recipes are convergent with 
the expected order and that the error in the entropy strategy is decreased by roughly a factor of 
two. 


The errors (Lı and Læ norms) for the four cases are shown in Fig. 6. Again, 
the second-order accuracy of the algorithm is recovered. Using the entropy strategy 
increases the numerical accuracy by around a factor of two and we suggest its use 
in the highly magnetised regime of BH magnetospheres. 


3.4 Magnetised equilibrium torus 

As a final validation of the code in the GRMHD regime, we perform the simulation 
of a magnetised equilibrium torus around a spinning BH. A generalisation of the 
steady-state solution of the standard hydrodynamical equilibrium torus with con- 
stant angular momentum [see, e.g., 43, 134, 135] to MHD equilibria with toroidal 
magnetic fields was proposed by [136]. This steady-state solution is important since 
it constitutes a rare case of a non-trivial analytic solution in GRMHD™!. 

For the initial setup of the equilibrium torus, we adopt a particular relationship 
w = w(p), where w = ph is the fluid enthalpy and & = (py), where Pm = Lpm, 
@ = Lw, Pm = b*/2 is the magnetic pressure, and L = gig9ie — 91196. From these 
relationships, thermal and magnetic pressures are described as 


p= Ku", (87) 
Pm = Kw". (88) 
The analytical solutions can be constructed from 


RPT Pm _ 
w Kue a as a a 


'8]Note that the discrepancy in v” appears less dramatic when viewed in terms of 


the four-velocity u”. 
91We thank Chris Fragile for providing subroutines for this test case. 
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Table 4 Parameters for the MHD equilibrium torus test 


Case lo Te W; We We Be 
A 2.8 4.62 -0.030 -0.103 1.0 0.1 


for the introduced total potential W, where W = In |u|. The centre of the torus 
is located at (re, 7/2). At this point, we parametrize the magnetic field strength in 
terms of the pressure ratio 


Be = Pg(Te,/2)/Pm(Te, 7/2) . (90) 
The gas pressure and magnetic pressure at the centre of the torus are given by 


K 1\ t 
Pe = we(Win We) (= 1 F T B ) , Pm. = Pef Be (91) 


From these, the constants K and K,,, for barotropic fluids are obtained. 
The magnetic field distribution is given by the distribution of magnetic pressure 
Pm. From the consideration of a purely toroidal magnetic field one obtains 


b? = /2pm/A, bt = Lb’, (92) 


where A = ggg + 2lgtg + l’ giz and £L := —ug/u is the specific angular momentum. 

We perform 2D simulations using logarithmic KS coordinates with h = 0 and 
Ro = 0. The simulation domain is 0 € [0,7], rks € [0.95rn, 50M], where ra is 
the (outer) event horizon radius of the BH. The BH has the dimensionless spin 
parameter a = 0.9. For simplicity, we set the two indices to the same value of 
k = ņ = 4/3 and also set the adiabatic index of the adopted ideal EOS to y = 4/3. 
The remaining parameters are listed in the Table 4. 

Initially, the velocity of the atmosphere outside of the torus is set to zero in 
the Eulerian frame, with density and gas pressure set to very small values of p = 
Pmin le pP = Pmin fo with pmin = 107° and pmin = 1077. It is important 
to note that the atmosphere is free to evolve and only densities and pressures are 
floored according to the initial state. In the simulations we use the HLL approximate 
Riemann solver, third order LIMO3 reconstruction, two-step time update, and a 
CFL number of 0.5. We impose outflow conditions on the inner and outer boundaries 
of the radial direction and reflecting boundary conditions in the 0 direction. As the 
magnetic field is purely toroidal, and will remain so during the time-evolution of 
this axisymmetric case, no particular V - B = 0 treatment is used. 

The top panels of Fig. 7 show the density distribution at the initial state and at t = 
200 M, as well as the plasma beta distribution at t = 200 M. The rotational period 
of the disk centre is tr = 68 M. The initial torus configuration is well maintained 
after several rotation period. For a qualitative view of the simulations, the 1D 
radial and azimuthal distributions of the density are shown in the lower two panels 
in Fig. 7 with different grid resolutions. All but the low resolution case are visually 
indistinguishable from the initial condition in the bottom-left panel, showing p — r 
with a linear scale. Since the atmosphere is evolved freely, small density waves 
propagate in the ambient medium of the torus, as seen in the p — 0 cut. This does 
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p, t=200,0M B, t =200.0M 


Figure 7 Top: qualitative view of the torus evolution at a resolution of Ny = Ng = 400. The 
spatial scale is given in units of M. Left: Initial rest-frame density distribution, center: density at 
t = 200 M, right: plasma 8 parameter at t = 200 M. Bottom: density slices through the torus at 
t = 100 M for constant 0 = 7/2 (left) and r = 5 (right). 


Error 


Figure 8 Error of the density p in the strongly-magnetised Komissarov torus, comparing the 
solution at t = 60 M with the exact solution. The second-order behaviour of the numerical 


scheme is well recovered. 
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not adversely affect the equilibrium solution in the bulk of the torus however. Error 
quantification (Lı and Lə) is provided in Fig. 8. The second-order properties of 
the numerical scheme are well recovered. 


3.5 Differences between FCT and GLM 

Having implemented two methods for divergence control, we took the opportunity 
to compare the results of simulations using both methods. We analysed three tests: 
a relativistic Orszag-Tang vortex, magnetised Michel accretion, and magnetised ac- 
cretion from a Fishbone-Moncrief torus. Although much less in-depth, this compar- 
ison is in the same spirit as those performed in previous works in non-relativistic 
MHD [100, 137, 138]. The well-known work by [100] compares seven divergence- 
control methods, including an early non-conservative divergence-cleaning method 
known as the eight-wave method [139], and three CT methods, finding that FCT 
is among the three most accurate schemes for the test problems studied. In [137], 
three divergence-cleaning schemes and one CT scheme were applied to the same 
test problem of supernova-induced MHD turbulence in the interstellar medium. It 
was found that the three divergence-cleaning methods studied suffer from, among 
other problems, spurious oscillations in the magnetic energy, which is attributed to 
the non-locality introduced by the loss of hyperbolicity in the equations. Finally, 
in [138], a non-staggered version of CT adapted to a moving mesh is compared to 
the divergence-cleaning Powell scheme [140], an improved version of the eight-wave 
method. They observe greater numerical stability and accuracy, and a better preser- 
vation of the magnetic field topology for the CT scheme. In their tests, the Powell 
scheme suffers from an artificial growth of the magnetic field. This is explained to 
be a result of the scheme being non-conservative. 


8.5.1 Orszag-Tang Vortex 

The Orszag-Tang vortex [141] is a common problem that can be used to test MHD 
codes for violations of V - B. The relativistic version presented here was performed 
in 2D using Cartesian coordinates in a 128 x 128-resolution domain of 27 x 27 length 
units with periodic boundary conditions, and evolved for 10 time units (c = 1). The 
equation of state was chosen to be that of an ideal fluid with 7 = 4/3 and the 
initial conditions were set to p = 1.0, p = 10.0, v” = —0.99siny, v” = 0.99 sin z, 
B®” = — siny and BY = sin2z . Snapshots of the evolution are shown in Fig. 9. 

As can be seen in Figs. 9 and 10, the general behaviour in both cases is quite 
similar qualitatively, with only slight differences at specific locations. For instance, 
when compared to GLM, FCT exhibited higher and sharper maxima for the mag- 
nitude of the magnetic field. In a similar fashion, some fine features in the Lorentz 
factor that can be seen in Fig. 9 for FCT appear to be smeared out when using 
GLM, giving a false impression of symmetry under 90° rotations, while the actual 
symmetry of the problem is under 180° rotations. This may be an evidence of FCT 
being less diffusive than GLM. 


3.5.2 Spherical accretion 
We tested the ability of both methods to preserve a stationary solution by evolving 
a magnetised Michel accretion in 2D, as shown in Fig. 11. We employed spherical 
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Figure 9 Relativistic Orszag-Tang vortex. Left column: small differences can be observed in this 
snapshot of the Lorentz factor at t = 5.0. Some features that appear when using CT are flattened 
when using GLM, possibly due to a greater diffusivity of the latter. Middle column: final snapshot 
of B*B;. Good agreement between the two methods can be seen, except at some extreme points. 
Right column: violation of V - B = 0. 


g10|Y -B/|B|| 


lo, 


logio|V -B/[B|| 


MKS coordinates (see Sect. 2.7.1), N, x Ng = 200 x 100 resolution, and a domain 
with r € [1.9M,10M] and 0 € [0,7]. The fluid obeyed an ideal equation of state 
with 4 = 5/3 and the sonic radius was located at re = 8, and the magnetic field 
was normalised so that the maximum magnetisation was 0 = 10°. We repeated 
the numerical experiment of section Sect. 3.3, now in 2D. As shown in Fig. 11, 
numerical artefacts start to become noticeable at these later times. For instance, 
with these extreme magnetisations, for GLM we observe spurious features near the 
poles at 0 = 0 and 0 = 7, as well as deviations in the velocity field near the outer 
boundary r = 10 M. The polar region is of special interest for jet simulations, where 
the divergence-control method must be robust enough to interplay with the axial 
boundary conditions. The bottom of Fig. 11 shows the profiles of several quanti- 
ties at 0 = 1/2. Both divergence-control methods produce an excellent agreement 
between the solution at different times in the equatorial region. The rightmost col- 
umn in the bottom of Fig. 11 shows the relative errors in the radial component of 
the magnetic field for each method. The errors for FCT are not only one order of 
magnitude lower than for GLM, but also behave differently, remaining at the same 
level near the more-magnetised inner region instead of growing as seen for GLM. 


3.5.8 Accreting Torus 

To compare both methods in a setting closer to our intended scientific applica- 
tions, we simulated accretion from a magnetised perturbed Fishbone-Moncrief torus 
around a Kerr BH with M = 1 and a = 0.9375. We employed modified spherical 
MKS coordinates as described in Sect. 2.7.1 and a domain where r € [1.29, 2500] 
and @ € [0,7] with a resolution of N, x Ng = 512 x 256, and evolved the system until 
t = 2000 M. At the radial boundaries, we imposed noinflow boundary conditions 
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— GM — GLM 
- FCT 


Figure 10 Relativistic Orszag-Tang vortex: cuts at y = 7/2 and t = 10.0 of the density p (left) 
and the magnitude of the magnetic field |B| (right). While for p there is in general a good 
agreement, FCT tends to produce higher maxima for the magnetic field. 


while at the boundaries with the polar axis we imposed symmetric boundary condi- 
tions for the scalar variables and the radial vector components and antisymmetric 
boundary conditions for the azimuthal and polar components. In the BHAC code, 
noinflow boundary conditions are implemented via continuous extrapolation of the 
primitive variables and by replacing the three-velocity with zero in case inflowing 
velocities are present in the solution. The fluid obeyed an ideal equation of state 
with 7 = 4/3. The inner edge of the torus was located at rin = 6.0 and the max- 
imum density was located at rmax = 12.0. The initial magnetic field configuration 
consisted of a single loop with Ag « (~/Pmax — Pcut) and zero for p < peut = 0.2. To 
simulate vacuum, the region outside the torus was filled with a tenuous atmosphere 
as is customarily done in these types of simulation. In this case, the prescription for 


=3/2 and Patm = Pmin r75/2, where Pmin = 1075 


the atmosphere was Patm = Pmin T 
and pmin = 1/3 x 1077. A qualitative difference can be seen even at early times 
of the simulation. The two upper panels of Fig. 12 show a snapshot of the simula- 
tion at t = 20 M using both GLM and FCT. For GLM some of the magnetic field 
has diffused out of the original torus, magnetising the atmosphere. This artefact is 
visible for GLM from almost the beginning of the simulation (t ~ 20 M), while for 
FCT it is minimal. Even though this particular artefact is not of crucial importance 
for the subsequent dynamics of the simulation, this points to a higher inclination of 
GLM to produce spurious magnetic field structures. At later times (bottom panels 
of Fig. 12), the most noticeable difference is the smaller amount of turbulent mag- 
netic structures and the bigger plasma magnetisation inside the funnel in FCT, as 
compared to GLM. This latter difference indicates that the choice of technique to 
control V - B may have an effect on the possibility of jet formation in GRMHD 
simulations, although this specific effect was not extensively studied. 

To summarise this small section on the comparison between both divergence- 
control techniques, we found from the three tests performed that FCT seems to be 
less diffusive than GLM, is able to preserve for a longer time a stationary solution, 
and seems to create less spurious structures in the magnetic field. However, it still 
has the inconvenient property that it is not possible to implement a cell-entered 
version of it whilst fully incorporating AMR. As mentioned previously, we are cur- 
rently working on a staggered implementation adapted to AMR, and this will be 
described in a separate work. 
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Figure 11 Top: logarithmic density and streamlines in 2D magnetised Michel accretion at times 

t =0M (left) and t = 100 M using GLM (centre) and FCT (right). The horizon is marked by the 
black line at r = 2. Bottom: profiles at 0 = 7/2 of, from left to right, radial 3-velocity, density 
and radial magnetic field at t = 0 M (blue circles) and t = 100 M (red line) using GLM (upper) 
and FCT (lower). The last column shows the relative difference between the magnetic field at 

t = 100 M and at the initial condition. 


4 Torus simulations 

4.1 Initial conditions 

We consider a hydrodynamic equilibrium torus threaded by a weak magnetic field 
loop. The particular equilibrium torus solution with constant angular momentum 
was first presented by [134] and [142] and is now a standard test for GRMHD 
simulations [see, e.g., 40, 50, 125, 135, 143]. To facilitate cross-comparison, we set 
the initial conditions in the torus close to those adopted by [67, 125]. Hence the 
spacetime is a Kerr BH with dimensionless spin parameter a = 0.9375. The inner 
radius of the torus is set to rin = 6 M and the density maximum is located at Tmax = 
12 M, where radial and azimuthal positions refer to Boyer-Lindquist coordinates. 
With these choices, the orbital period of the torus at the density maximum becomes 
T = 247M. We adopt an ideal gas EOS with an adiabatic index of 4 = 4/3. A weak 
single magnetic field loop defined by the vector potential 


Ag x max(p/Pmax — 0.2,0), (93) 


is added to the stationary solution. The field strength is set such that 2pmax/b7 4. = 


max 
100, where global maxima of pressure pmax and magnetic field strength b?,,, do 
not necessarily coincide. In order to excite the MRI inside the torus, the thermal 
pressure is perturbed by 4% white noise. 

As with any fluid code, vacuum regions must be avoided and hence we apply 


floor values for the rest-mass density (pg = 10~°r~*/?) and the gas pressure (pg = 
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Figure 12 Magnetised torus: plasma 6 at t = 20 M (top) and density and magnetic field lines 
t = 2000 M (bottom) using GLM (left) and FCT (right). 


1/3 x 1077 r75/2). In practice, for all cells which satisfy p < pa we set p = pa, in 
addition if p < pg, we set p = pa. 

The simulations are performed using horizon penetrating logarithmic KS coor- 
dinates (corresponding to our set of modified KS coordinates with h = 0 and 
Ro = 0). In the 2D cases, the simulation domain covers rgs € [0.96rp, 2500 M] 
and 0 € [0,7], where rp ~ 1.35 M. In the 3D cases, we slightly excise the axial re- 
gion 6 € [0.027, 0.987] and adopt ¢ € [0,27]. We set the boundary conditions in the 
horizon and at r = 2500 M to zero gradient in primitive variables. The 6-boundary 
is handled as follows: when the domain extends all the way to the poles (as in our 
2D cases), we adopt “hard” boundary conditions, thus setting the flux through the 
pole manually to zero. For the excised cone in the 3D cases, we use reflecting “soft” 
boundary conditions on primitive variables. 

The time-update is performed with a two-step predictor corrector based on the 
TVDLF fluxes and PPM reconstruction. Furthermore, we set the CFL number to 
0.4 and use the FCT algorithm. Typically, the runs are stopped after an evolution 
for t = 5000 M, ensuring that no signal from the outflow boundaries can disturb the 
inner regions. To check convergence, we adopt the following resolutions: N, x Ng € 
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{256 x 128,512 x 256, 1024 x 512} in the 2D cases and N, x Ng x Ng € {128 x 
64 x 64,192 x 96 x 96,256 x 128 x 128,384 x 192 x 192} in the 3D runs. In the 
following, the runs are identified via their resolution in 9-direction. For the purpose 
of validation, we ran the 2D cases also with the HARM3D code [77].!!° 

To facilitate a quantitative comparison, we report radial profiles of disk-averaged 
quantities similar to [67, 125, 144]. For a quantity q(r,0,¢,t), the shell average is 
defined as 


27 


lo(r,t)) = I2 Semin 4". 9, & v= 9 do a 


m pOmax 
k Omin V =g do dð 


(94) 


which is then further averaged over a given time interval to yield (q(r)) (note that we 
omit the weighting with the density as done by [67, 125]). The limits Onin = 7/3, 
Omax = 27/3 ensure that atmosphere material is not taken into account in the 
averaging. The time-evolution is monitored with the accretion rate M and the 
magnetic flux threading the horizon ¢g 


27 T 
M := r /=gd0 do, 95 
| [my g dô do (95) 
1 27 T E 
bs=3 f f 1B"V=5 ddo, (96) 


where both quantities are evaluated at the outer horizon rp. 


4.2 2D results 

Figure 13 illustrates the qualitative time evolution of the torus by means of the rest- 
frame density p, plasma-3 and the magnetisation o = b?/p. After t ~ 300 M, the 
MRI-driven turbulence leads to accretion onto the central BH. The accretion rate 
and magnetic flux threading the BH then quickly saturate into a quasi-stationary 
state (see also Fig. 14). The accreted magnetic flux fills the polar regions and gives 
rise to a strongly magnetised funnel with densities and pressures dropping to their 
floor values. For the adopted floor values we hence obtain values of plasma-3 as low 
as 107% and magnetisations peaking at o ~ 10? in the inner BH magnetosphere. 
These extreme values pose a stringent test for the robustness of the code and, 
consequently, the funnel region must be handled with the auxiliary inversion based 
on the entropy switch (see Sect. 2.10.2). 


4.2.1 Comparison to HARM3D 


For validation purposes we simulated the same initial conditions with the HARM3D code. 


Wherever possible, we have made identical choices for the algorithm used in both 
codes, that is: PPM reconstruction, TVDLF Riemann solver and a two step time 
update. It is important to note that the outer radial boundary differs in both codes: 
while the HARM3D setup implements outflow boundary conditions at r = 50M, in 
the BHAC runs the domain and radial grid is doubled in the logarithmic Kerr-Schild 
coordinates, yielding identical resolution in the region of interest. This ensures that 


10 The results were kindly provided by Monika Moscibrodzka. 
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Figure 13 Evolution of the 2D magnetised torus with resolution 1024 x 512 for times 

t/M € {300, 1000, 2000}. We show logarithmic rest-frame density (top), logarithmic plasma 8 
(middle) and the logarithm of the magnetisation parameter o = b? /p (bottom). Magnetic field 
lines are traced out in the first panel using black contour lines. One can clearly make out the 
development of the MRI and evacuation of a strongly magnetised funnel reaching values of 

B < 10-5 and ø ~ 10°. 
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no boundary effects compromise the BHAC simulation. Next to the boundary condi- 
tions, also the initial random perturbation varies in both codes which can amount 
to a slightly different dynamical evolution. 

After verifying good agreement in the qualitative evolution, we quantify with both 
codes M and ¢g according to equations (95) and (96). The results are shown in 
Fig. 14. Onset-time of accretion, magnitude and overall behaviour are in excellent 
agreement, despite the chaotic nature of the turbulent flow. We also find the same 
trend with respect to the resolution-dependence of the results: upon doubling the 
resolution, the accretion rate (M ), averaged over t € [1000, 2000], increases signif- 
icantly by a factor of 1.908 and 1.843 for BHAC and HARM , respectively. For (dp), 
the factors are 1.437 and 1.484. At a given resolution, the values for (M) and (dg) 
agree between the two codes within their standard deviations. Furthermore, we have 
verified that these same resolution variations are within the run-to-run deviations 
due to a different random number seed for the initial perturbation. 


— BHAC, 256 — HARM3D, 256 


— BHAC, 512  — HARM3D, 512 
0 500 1000 1500 2000 
t M] 


Figure 14 Accretion rates and horizon-penetrating magnetic flux in the 2D validation runs. We 
show two resolutions with each code: BHAC (blue, red) and HARM3D (dark blue, orange). Despite 
the chaotic nature of the turbulent accretion both codes show very good qualitative and 
quantitative agreement. 


Further validation is provided in Fig. 15 where disk-averaged profiles for the two 
highest resolution 2D runs are shown according to equation (94). The quantities of 
interest are the rest-frame density p, the dimensionless temperature © := p/pc?, 
the magnitude of the fluid-frame magnetic field |B| = Vb, thermal and mag- 
netic pressures Peas, Pmag and the plasma-G. Again we set the averaging time 
t € [1000, 2000] M with both codes. The agreement can be considered as very good, 
that is apart from a slightly higher magnetisation in HARM for r € [20,30], the 
differences of which are well within the lo standard deviation over the averaging 
time. Small systematic departures at the outer edge of the HARM domain are likely 
attributable to boundary effects. 


4.3 3D results 
We now turn to the 3D runs performed with BHAC. The qualitative evolution of the 
high resolution run is illustrated in Fig. 16 showing rest-frame density and b? on the 
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Figure 15 Disk-averaged quantities in the 2D validation runs. The blue curves are obtained with 
BHAC and the red curves with HARM3D in a two-dimensional setting. The shaded regions mark the 
lo standard deviation of the spatially-averaged snapshots (omitted for the highly fluctuating (3)). 
Apart from a slightly higher magnetisation in HARM for r € [20, 30], we find excellent agreement 
between both codes. 


two slices z = 0 and y = 0. Overall, the evolution progresses in a similar manner 
to the 2D cases: MRI-driven accretion starts at t ~ 300 M and enters saturation 
at around t ~ 1000 M. Similar values for the magnetisation in the funnel region 
are also obtained. However, since the MRI cannot be sustained in axisymmetry as 
poloidal field cannot be re-generated via the ideal MHD induction equation [145], 
we expect to see qualitative differences between the 2D and 3D cases at late times. 

Four different numerical resolutions were run which allows a first convergence anal- 
ysis of the magnetised torus accretion scenario. Based on the convergence study, we 
can estimate which numerical resolutions are required for meaningful observational 
predictions derived from GRMHD simulations of this type. 

Since we attempt to solve the set of dissipation-free ideal MHD equations, conver- 
gence in the strict sense cannot be achieved in the presence of a turbulent cascade 
[see also the discussion in 146, 147]. Instead, given sufficient scale separation, 
one might hope to find convergence in quantities of interest like the disk averages 
and accretion rates. The convergence of various indicators in similar GRMHD torus 
simulations was addressed for example by [67]. The authors found signs for con- 
vergence in most quantifications when adopting a resolution of 192 x 192 x 128, 
however no convergence was found in the correlation length of the magnetic field. 
Hence the question as to whether GRMHD torus simulations can be converged with 
the available computational power is still an open one. 

From Figs. 17 and 18, it is clear that the resolution of the Ng = 64 run is insuf- 
ficient: a peculiar mini-torus is apparent in the disk-averaged density which dimin- 


n Even when the dissipation length is well resolved, high-Reynolds number flows 
show indications for positive Lyapunov exponents and thus non-convergent chaotic 
behaviour [see, e.g., 148]. 
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Figure 16 Fluid-frame density (top) and log,, b? (bottom) for t = 3000 M on the y = 0 plane 
(left) and the z = 0 plane (right) in the 3D magnetised torus run with resolution 384 x 192 x 192. 
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ishes with increasing resolution. Also the onset-time of accretion and the saturation 
values differ significantly between the Ng = 64 run and its high-resolution counter- 
parts. These differences diminish between the high-resolution runs and we can see 
signs of convergence in the accretion rate: increasing resolution from Ng = 128 to 
No = 192 appears to not have a strong effect on M. Also the evolution of ¢g agrees 
quite well between Ng = 128 and Ng = 192. Hence the systematic resolution depen- 
dence of M and ¢g in the (even higher resolution) 2D simulations appears to be an 
artefact of the axisymmetry. It is also noteworthy that the variability amplitude of 
the accretion rate is reduced in the 3D cases. It appears that the superposition of 
uncorrelated accretion events distributed over the ¢-coordinate tends to smear out 
the sharp variability that results in the axisymmetric case. 


10° 


= 101 


(0) 1000 2000 3000 4000 5000 
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Figure 17 Accretion rates and horizon-penetrating magnetic flux in the 3D runs for varying 
numerical resolution. We show results from four different resolutions labeled according to the 
number of cells in 6-direction. 


Although the simulations were run until t = 5000 M, in order to enable com- 
parison with the 2D simulations, we deliberately set the averaging time to t € 
[1000 1/7, 2000 M]. Figure 18 shows that as the resolution is increased, the disk- 
averaged 3D data approaches the much higher resolution 2D results shown in Fig. 15, 
indicating that the dynamics are dominated by the axisymmetric MRI modes at 
early times. When the resolution is increased from Ng = 128 to Nọ = 192, the 
disk-averaged profiles generally agree within their standard deviations, although 
we observe a continuing trend towards higher gas pressures and magnetic pres- 
sures in the outer regions r € [30M, 50M]. The overall computational cost quickly 
becomes significant: for the Ng = 128 simulation we spent 100 K CPU hours on 
the Iboga cluster equipped with Intel(R) Xeon(R) E5-2640 v4 processors. As the 
runtime scales with resolution according to Nj, doubling resolution would cost a 
considerable 1.6 M CPU hours. 


4.4 Effect of AMR 
In order to investigate the effect of the AMR treatment, we have performed a 2D 
AMR-GRMHD simulation of the torus setup. It is clear that whether a simulation 
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Figure 18 Disk-averaged quantities in the 3D runs for varying numerical resolution. The shaded 
regions mark the lo standard deviation of the spatially-averaged snapshots as in Fig. 15. 


can benefit from adaptive mesh refinement is very much dependent on the physical 
scenario under investigation. For example, in the hydrodynamic simulations of re- 
coiling BHs due to [37], refinement on the spiral shock was demonstrated to yield 
significant speedups at a comparable quality of solution. This is understandable 
as the numerical error is dominated by the shock hypersurface. In the turbulent 
accretion problem, whether automated mesh refinement yields any benefits is not 
clear. 

The initial conditions for this test are the same as those used in Sect. 4.1. However, 
due to the limitation of current AMR treatment, we resort to the GLM divergence 
cleaning method. Three refinement levels are used and refinement is triggered by 
the error estimator due to [127] with the tolerance set to e, = 0.1 (see Sect. 2.11). 
The numerical resolution in the base level is set to N, x Ng = 128 x 128. To test 
the validity and efficiency, we also perform the same simulation in a uniform grid 
with resolution of N, x Ng = 512 x 512 which corresponds to the resolution on the 
highest AMR level. 

Figure 19 shows the densities at t = 2000 M as well as the time-averaged density 
and plasma beta for the AMR and uniform cases. The averaged quantities are cal- 
culated in the time interval of t € [1000 M, 2000 M]. The overall behaviour is quite 
similar in both cases. Naturally, differences are seen in the turbulent structure in the 
torus and wind region for a single snapshot. However, in terms of averaged quan- 
tities, the difference becomes marginal. In order to better quantify the difference 
between the AMR and uniform runs, the mass accretion rate and horizon penetrat- 
ing magnetic flux are shown in Fig. 20. These quantities exhibit a similar behaviour 
in both cases. In particular, the difference between the AMR run and the uniform 
run is smaller than the one from different resolution uniform runs and compatible 
with the run-to-run variation due to a different random number seed (cf. Sect. 4.2). 


This is unsurprising since the error estimator triggers refinement of the innermost 
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Figure 19 2D logarithmic density at t = 2000 M (left), averaged density (middle), and averaged 
plasma beta (right) of the 2D magnetised torus with three-levels AMR (top panels) and uniform 
resolution 512 x 512 (bottom panels). Magnetic field lines are traced out in the middle panels 
using black contour lines. The averaged quantities are calculated in the time interval 

t € [1000 M, 2000 M]. AMR blocks containing 16? cells are indicated in the upper left panel. 
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torus region to the highest level of AMR during most of the simulation time. The 
development of small scale turbulence by the MRI is clearly captured and it leads 
to similar mass accretion onto the BH. 


0 500 1000 1500 2000 
tM] 


Figure 20 Accretion rates and horizon penetrating magnetic flux of the 2D magnetised torus with 
three levels of AMR (black) and uniform resolution 512 x 512 (red). 


Table 5 CPU hours (CPUH) spent by the simulations of the 2D magnetised torus at uniform 
resolution and fraction of that time spent by the equivalent AMR runs up to t = 2000 M. 


Grid size CPU time Equiv. AMR 


(Nr x No) uniform time fraction 
[CPUH] [es = 0.1] 
512 x 512 674.0 0.643 


One of the important merits of using AMR is the possibility to resolve small and 
large scale dynamics simultaneously with lower computational cost than uniform 
grids. Figure 21 shows the large scale structure of the averaged magnetisation after 
10000 M of simulation time. The averaged quantities are calculated in the time 
interval ¢ € [6000 M, 10000 M]. In order to allow the large-scale magnetic field 
structure to settle down, we average over a later simulation time compared to the 
previous non-AMR cases. From the figure the collimation angle and magnetisation 
of the highly magnetised funnel in the AMR case are slightly wider than those in 
uniform case but the large-scale global structure is very similar in both cases. 

A comparison of the computational time for a uniform resolution with 512? and 
the equivalent AMR run (three-level AMR) is shown in Table 5. It is encouraging 
that even in the naive three-level AMR simulation we obtain qualitatively similar 
results comparable to the high resolution uniform run, but with having spent only 
64% of the computational time of the uniform run.!!?) Figure 22 shows the evolution 


121 Since we use the same Courant limited timestep for all grid-levels, the speedup 
is entirely due to saving in computational cells. The additional speedup that would 
be gained from [149]-type hierarchical timesteps can be estimated from the level 
population of the simulation: the expected additional gain is only ~ 8% for this 
setup. 
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Figure 21 2D logarithmic averaged magnetisation of the magnetised torus with three levels of 
AMR (left) and uniform resolution 512 x 512 (right). Magnetic field lines are traced out by white 
contour-lines. The averaged quantities are calculated in the time interval of 

t € [6000 M, 10000 M]. 
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Figure 22 Number of cells as a function of time for the AMR simulation. The dotted lines show 
the resolution of uniform grids equivalent to each of the three AMR levels. 
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of the total number of cells during the simulations of AMR cases. Initially less than 
216 cells are used even when we use three AMR levels, which is a similar number 
of cells as the uniform grid case with 256 x 256. When the simulation starts, the 
total cell number increases rapidly due to development of turbulence in the torus 
which is triggering higher refinement. We note that the total number of cells is 
still half of the total number of cells in the corresponding high-resolution uniform 
grid simulation (512 x 512), thus resulting in a direct reduction of computational 
cost. With increasing dynamic range, we expect the advantages of AMR to increase 
significantly, rendering it a useful tool for simulations involving structures spanning 
multiple scales. We leave a more detailed discussion on the effect of the AMR 
refinement strategy and various divergence-control methods to a future paper. 


5 Radiation post-processing 

In order to compute synthetic observable images of the BH shadow and surrounding 
accretion flow it is necessary to perform general-relativistic ray-tracing and GRRT 
post-processing [see, e.g., 73, 150-156]. In this article the GRRT code BHOSS (Black 
Hole Observations in Stationary Spacetimes) [78] is used to perform these calcula- 
tions. From BHAC, GRMHD simulation data are produced which are subsequently 
used as input for BHOSS. Although BHAC has full AMR capabilities, for the GRRT 
it is most expedient to output GRMHD data that has been re-gridded to a uniform 
grid. 

Since these calculations are performed in post-processing, the effects of radia- 
tion forces acting on the plasma during its magnetohydrodynamical evolution are 
not included. Additionally, the fast-light approximation has also been adopted in 
this study, i.e., it is assumed that the light-crossing timescale is shorter than the 
dynamical timescale of the GRMHD simulation and the dynamical evolution of 
the GRMHD simulation as light rays propagate through it is not considered. Such 
calculations are considered in an upcoming article [78]. 

Several different coordinate representations of the Kerr metric are implemented in 
BHOSS, including Boyer-Lindquist (BL), Logarithmic BL, Cartesian BL, Kerr-Schild 
(KS), Logarithmic KS, Modified KS and Cartesian KS. All GRMHD simulation data 
used in this study are specified in Logarithmic KS coordinates. Although BHOSS can 
switch between all coordinate systems on the fly, it is most straightforward to per- 
form the GRRT calculations in the same coordinate system as the GRMHD data, 
only adaptively switching to e.g., Cartesian KS when near the polar region. This 
avoids the need to transform between coordinate systems at every point along every 
ray in the GRMHD data interpolation, saving computational time. 


5.1 Radiative transfer equation 

Electromagnetic radiation is described by null geodesics of the background space- 
time (in this case Kerr), and these are calculated in BHOSS using a Runge-Kutta- 
Fehlberg integrator with fourth order adaptive step sizing and 5th order error con- 
trol. Any spacetime metric may be considered in BHOSS, as long as the contravariant 
or covariant metric tensor components are specified, even if they are only tabulated 
on a grid. For the calculations presented in this article the Kerr spacetime is written 
algebraically and in closed-form. 
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The observer is calculated by constructing a local orthonormal tetrad using trial 
basis vectors. These basis vectors are then orthonormalized using the metric tensor 
through a modified Gram-Schmidt procedure. The initial conditions of each ray for 
the coordinate system under consideration are then calculated and the geodesics 
are integrated backwards in time from the observer, until they either: (i) escape to 
infinity (exit the computational domain), (ii) are captured by the BH, or (iii) are 
effectively absorbed by the accretion flow. 

In order to perform these calculations the GRRT equation is integrated in parallel 
with the geodesic equations of motion of each ray. Written in covariant form, the 
(unpolarized) GRRT equation, in the absence of scattering, may be written [152] 
as 


dI iy 
C (-avo7 + ža) | (97) 


where Z := I,,/v® is the Lorentz-invariant intensity, I, is the specific intensity, v is 
the frequency of radiation, a,,9 is the specific absorption coefficient and j,,9 is the 
specific emission coefficient. The subscript “v” denotes evaluation of a quantity at 
a specific frequency, v, and a subscript “0” denotes evaluation of a quantity in the 
local fluid rest frame. The terms k, and u” are the photon 4-momentum and the 
fluid 4-velocity of the emitting medium, respectively. The former is calculated from 
the geodesic integration and the latter is determined from the GRMHD simulation 
data. The affine parameter is denoted by A. 
By introducing the optical depth along the ray 


A 
BOS f wia, (98) 


0 


together with the Lorentz-invariant emission coefficient (n = jy v?) and Lorentz- 
invariant absorption coefficient (x = vap), the GRRT equation (97) may be rewrit- 
ten as 


L 
Wy ea (99) 
dt, x 


Following [152], equation (99) may be reduced to two differential equations 


dt, 


See an = TIE 1 
TI Qv, 0 ( 00) 
dI — Jv,0 
ID = ve exp ( Tv) ; (101) 
where 
kau% 
V _ (kau) ons (102) 


y= = TkpuP) ; 


is the relative energy shift between the observer (“obs”) and the emitting fluid 
element. Integrating the GRRT equation in terms of the optical depth in the manner 
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presented provides two major advantages. Firstly, the calculation of the geodesic 
and of the radiative transfer equation may be performed simultaneously, rather than 
having to calculate the entire geodesic, store it in memory, and then perform the 
radiative transfer afterwards. Secondly, by integrating in terms of the optical depth 
we may specify a threshold value (typically of order unity) whereby the geodesic 
integration is terminated when encountering optically thick media exceeding this 
threshold. The combination of these two methods saves significant computational 
time and expense. 


5.2 BHOSS-simulated emission from Sgr A* 

Having in mind the upcoming radio observations of the BH candidate Sgr A* at the 
Galactic Centre, the following discussion presents synthetic images of Sgr A*. The 
GRMHD simulations evolve a single fluid (of ions) and are scale-free in length and 
mass. Consequently a scaling must be applied before performing GRRT calculations. 
Within BHOSS this means specifying the BH mass, which sets the length and time 
scales, and specifying either the mass accretion rate or an electron density scale, 
which scales the gas density, temperature and magnetic field strength to that of a 
radiating electron. 

Since the GRMHD simulation is of a single fluid, it is necessary to adopt a prescrip- 
tion for the local electron temperature and rest-mass density. Several such prescrip- 
tions exist, some which scale using the mass accretion rate [see, e.g., 71, 157, 158], 
scale using density to determine the electron number density and physical accretion 
rate [see, e.g., 73, 159], and some by employing a time-dependent smoothing model 
of the mass accretion rate [see, e.g., 67]. 

The dimensionless proton temperature, Op, is defined as 

és." (103) 


2° 
Mpc 


where kg is the Boltzmann constant, T, is the geometrical (i.e., in physical units) 
proton temperature and mp is the proton mass. This is then calculated from the 
GRMHD simulation density (p) and pressure (p) as 


o= +, (104) 


where the fact that the equation of state is ideal and that 7 = 4/3 has been assumed. 
The magnetic field strength in geometrical units, Bgeo, is readily obtained from the 
code magnetic field strength B = ,/b,,b! as 


1/2 
Beco =€ (=) B. (105) 
p 


What remains is to specify T, (or Oe := kpT./mec”) and pgeo. For simplicity we 
adopt the prescription of [71], wherein Tp/Te is assumed to be a fixed ratio. Whilst 
such an approximation is rather crude, to zeroth order the protons and electrons 
may be assumed to be coupled in this way. To scale the electron number density 
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we adopt the method of [73], assuming a density scale typically of order 107 cm~. 
A somewhat more sophisticated approach is to employ a thresholding of the fluid 
plasma beta where, when the local plasma beta exceeds some threshold the elec- 
trons and protons are coupled as previously mentioned (disk region), but when not 
exceeded (typically in the funnel region) the electron temperature is assumed to be 
constant [73, 76, 157]. Since plasma beta is found to decrease with resolution [67] 
and in this paper we seek only to demonstrate the convergence of our simulated 
shadow images obtained from the GRMHD data in regions where the density is 
non-negligible, we adopt the former model. 

For the plasma emissivity we use the approximate formula for thermal magneto- 
bremsstrahlung [160], which is given by 


. V2ne? Us 1/2 11/12 y1/6 2 1/3 
z T AC (x +21/2 X ) exp (-X ie (106) 


where e is the electron charge, ne the electron number density, and 


V 
= — 1 
X 5 (107) 
uv = ( 6 ) Božsinv, (108) 
ITMeC 


and v is the pitch angle of the photon with respect to the magnetic field. The 
absorption coefficient is readily obtained from Kirchoff’s law. 

Each image is generated using a uniform grid of 1000 x 1000 rays, sampling 60 
uniformly logarithmically spaced frequency bins between 109 Hz and 10!° Hz. All 
panels depict the observed image as seen at an observational frequency of 230 GHz, 
i.e. the frequency at which the EHT will image Sgr A*. This resolution is chosen 
because the integrated flux over the entire ray-traced image is convergent: doubling 
the resolution from 500 x 500 to 1000 x 1000 yields an increase of 0.17%, and from 
1000 x 1000 to 2000 x 2000 an increase of 0.09%. 

In practical GRRT calculations only simulation data which has already reached 
a quasi-steady state, typically t > 2000 M, is used. In this study we focus on the 
observational appearance of the accretion flow and BH shadow image. The detailed 
discussion of the spectrum, variability and plasma models warrants a separate study. 


5.3 Comparison of images 
A natural and important question arises from GRRT calculations of BH shadows: 
do ray-traced images of GRMHD simulation data converge as the resolution of the 
GRMHD simulation is increased? The existence of an optimal resolution, beyond 
which differences in images are small, implies that one can save additional com- 
putational time and expense by running the simulation at this optimal resolution. 
It would also imply that the GRMHD data satisfactorily capture the small-scale 
structure, turbulence and variations of the accretion flow. As such, it is informa- 
tive to investigate the convergence of BH shadow images obtained from GRMHD 
simulation data of differing resolutions, both quantitatively and qualitatively. 

To address this question we first generate a series of four snapshot images at 
t = 2500 M of the the accretion flow and BH shadow from GRMHD simulation 
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Figure 23 Snapshot images of 3D GRMHD simulation data with parameters chosen to mimic the 


emission from Sgr A*. The resolution of the simulation data is indicated in the bottom-right 
corner of each panel and discussed in the text. 
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data. The resolution of these data are 2M x N x N in (r,@, ¢), i.e., twice as much 
resolution in the radial direction compared to the zenith and azimuthal directions. 
The images depicted in Fig. 23 correspond to N = 64, 96, 128 and 192 respectively. 
Here, the proton to electron temperature ratio was chosen as T/T. = 3 (similar 
to [71, 157]), the electron number density scaling as 5 x 10’ cm~3, the BH mass is 
set to 4.5 x 10° Mo, the source distance is 8.4 x 10° pc, the dimensionless BH spin 
parameter 0.9375 and the observer inclination angle with respect to the BH spin 
axis is 60°. 

A direct consequence of increasing the resolution of the GRMHD data is resolv- 
ing the fine-scale turbulent structure of the accretion flow. The characteristic dark 
shadow delineating the BH shadow can be clearly seen in all images. As the resolu- 
tion of the GRMHD data is increased, the images become less diffuse. It is difficult 
with the naked eye to draw firm physical conclusions, and so in the following we 
perform a quantitative pixel-by-pixel analysis. 

With these snapshot images we may perform a quantitative measure of the 
difference between any two images through introducing the (normalised) cross- 
correlation. For two given two-dimensional arrays f(x,y) and g(x,y) (i.e., 2D im- 
ages), a measure of similarity or difference may be calculated through the cross- 
correlation C, where C € [—1,1]. The normalised cross-correlation is defined as 


1 


C=C, = 
Nofog 


tJ 


NC {LF (@,9) — usl lola, y) — nal} ; (109) 


where uf, of and ug, Og correspond to the mean and standard deviation of f and g 
respectively, and N is equal to the size of either f or g. In the examples considered 
here the images are all of equal size and dimension, so N = Ny = Ng. Equation 
(109) may be interpreted as the inner product between two data arrays, with the 
value of C expressing the degree to which the data are aligned with respect to each 
other. When C = 1 the data are identical, save for a multiplicative constant, when 
C = 0 the data are completely uncorrelated, and when C < 0 the data are negatively 
correlated. 

Each image pixel has an intensity value represented as a single greyscale value 
between zero and one. Given the relative intensity data of two different images, 
Equation (109) is then employed to determine the normalised cross-correlation be- 
tween the two images. This procedure applied to the panels in Fig. (23) yields the 
following symmetric matrix of cross-correlation values between the images 


1 0.839495 0.809205 0.856958 
fies | 1 0.867578 0.917560 (110) 
cae S = 1 0.948544 |` 

= = = 1 


Indices i and j, where (i, j) = (1,4), denote the images being cross-correlated. 
The rightmost column of equation (110) denotes the cross-correlation values, C;,4, 
in descending order between images, i.e., the cross-correlation of image 4 with im- 
ages 1, 2, 3 and 4 respectively. Since Ci+1,4 > Cj,4 it is clear that the similarity 
between images increases as the resolution of the GRMHD simulation is increased. 
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Similarly, for image 3 it is found that Ci+1,3 > C;,3. Finally, it also follows that 
C34 > C2,3 > C1,2, i.e., the correlation between successive pairs of images increases 
with increasing resolution, demonstrating the convergence of the GRMHD simula- 
tions with increasing grid resolution. Whilst the lowest resolution of 128 x 64 x 64 is 
certainly insufficient, both difference images and cross-correlation measures indicate 


that a resolution of 256 x 128 x 128 is sufficient and represents a good compromise. 


Figure 24 Matrix of image differences D; j of the four panels in Fig. 23. Upper diagonal panels 
are greyscale differences. Lower diagonal panels are identical to corresponding upper diagonal 
panels but with differences illustrated with RGB pixel values. Black panels correspond to D;,;, 
i.e., trivially the difference between an image and itself. 


6 Conclusions and outlook 

We have described the capabilities of BHAC, a new multidimensional general- 
relativistic magnetohydrodynamics code developed to perform hydrodynamical and 
MHD simulations of accretion flows onto compact objects in arbitrary stationary 
spacetimes exploiting the numerous advantages of AMR techniques. The code has 
been tested with several one-, two- and three- dimensional scenarios in special- 
relativistic and general-relativistic MHD regimes. For validation, GRMHD simu- 
lations of MRI unstable tori have been compared with another well-known and 
tested GRMHD code, the HARM3D code. BHAC shows very good agreement with the 
HARM3D results, both qualitatively and quantitatively. As a first demonstration of 
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the AMR capabilities in multi-scale simulations, we performed the magnetized- 
torus accretion test with and without AMR. Despite the latter intrinsically implies 
an overhead of ~ 10%, the AMR runtime amounted to 65% of that relative to the 
uniform grid, simply as a result of the more economical use of grid cells in the 
block based AMR. At the same time, the AMR results agree very well with the 
more expensive uniform-grid results. With increasing dynamic range, we expect the 
advantages of AMR to increase even more significantly, rendering it a useful tool 
for simulations involving structures of multiple physical scales. 

Currently, two methods controlling the divergence of the magnetic field are avail- 
able in BHAC and we compared them in three test problems. Although solutions 
obtained with the cell-centered flux-interpolated constrained transport (FCT) algo- 
rithm and the divergence cleaning scheme (GLM) yield the same (correct) physical 
behaviour in the case of weak magnetic fields, FCT performs considerably better in 
the presence of strong magnetic fields. In particular, FCT is less diffusive than GLM, 
is able to preserve a stationary solution, and it creates less spurious structures in 
the magnetic field. For example, the use of GLM in the case of accretion scenarios 
with strong magnetic fields leads to worrisome artefacts in the highly magnetised 
funnel region. The development of a constrained transport scheme compatible with 
AMR is ongoing and will be presented in a separate work [161]. 

The EHTC and its European contribution, the BlackHoleCam project [68], aim 
at obtaining horizon-scale radio images of the BH candidate at the Galactic Cen- 
ter. In anticipation of these results, we have used the 3D GRMHD simulations 
as input for GRRT calculations with the newly developed BHOSS code [78]. We 
found that the intensity maps resulting from different resolution GRMHD simu- 
lations agree very well, even when comparing snapshot data that was not time 
averaged. In particular, the normalised cross-correlation between images achieves 
up to 94.8% similarity between the two highest resolution runs. Furthermore, the 
agreement between two images converges as the resolution of the GRMHD simula- 
tion is increased. Based on this comparison, we find that moderate grid resolutions 
of 256 x 128 x 128 (corresponding to physical resolutions of Arxg x AOxg x Adxg = 
0.04M x 0.024rad x 0.05rad at the horizon) yield sufficiently converged intensity 
maps. Given the large and likely degenerate parameter space and the uncertainty 
in modelling of the electron distribution, this result is encouraging, as it demon- 
strates that the predicted synthetic image is quite robust against the ever-present 
time variability, but also against the impact that the grid resolution of the GRMHD 
simulations might have. In addition, independent information on the spatial orien- 
tation and magnitude of the spin, such as the one that could be deduced from the 
dynamics of a pulsar near Sgr A* [162], would greatly reduce the space of degener- 
ate solutions and further increase the robustness of the predictions that BHAC will 
provide in terms of synthetic images. 

Finally, we have demonstrated the excellent flexibility of BHAC with a variety of 
different astrophysical scenarios that are ongoing and will be published shortly. 
These include: oscillating hydrodynamical equilibrium tori for the modelling of 
quasi-periodic oscillations [163], episodic jet formation [164] and magnetised tori 
orbiting non-rotating dilaton BHs [105]. 
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Appendix A: Evolution of the scalar @ 


To obtain the evolution equation for ¢ in the augmented Faraday’s law, we project 
(42) onto the Eulerian observer by contracting with —n,, as 


-VCF ny, — bn’) = =r — (CF HY — dg!” )V yny (111) 

=> V B” + Von’ = —Kb — FH” — bg°”)Viny (112) 

=> (-9) Pah PaB] + Vion” = -ro — CFH — bg”) Ving , (113) 
where we used BY = —n, *F#”. Using the definition of extrinsic curvature Ky, := 


=à Vany, we can write [Eq. (7.62) in 40] 

Vinu = Nyau — Kw, (114) 
where we used the “acceleration” of the Eulerian observer a, := mV nN, which sat- 
isfies na, = 0. With the identity a; = a™ t;a [165] and exploiting the symmetries 
of *F”” and Ky,, is straightforward to show that 

(Pala) aB] + oF’ Vn, = (7) Palla) aB] — Bia. (115) 


Hence it follows that 


d(T) + iL 1(aB’ — 96") = — aP rg + ay? pg” Vn, 


116 
+7! Baa. l ) 


Using again Eq. (114), the source term S := \/ya¢g"”V_n, can be rewritten as 


S= =,/ yaa” Kuv, (117) 


where the first term drops out due to the orthogonality na, = 0. For a symmetric 


tensor S#”, we have 
= of ney tee: 
aS” Kuv = aS" Kij = S50;B? + zS O" Onis . (118) 


This follows from the relation T9; = —Kjja~' where TP), are elements of the 4- 
Christoffel symbols [see e.g., (B.9) of 85]. Thus 


S =~ V74(018" + 577 Ours) (119) 


= ~V7bdi6' — VIPE Y HOR (120) 
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Appendix B: Modified Faraday’s law 


The augmented Faraday’s law follows from the j-component of (42) as 


VCE” — 697”) = kopi /a (121) 
=(—g)~/? {0:[,/7(—B*)] + 0; [Vy (VIB? — V'B?)]} + gO) (—¢) = kop! /a 
(122) 


=ô; (VIBİ) + 0; [V7 (V'B? — VI B’)]| + Vyag?* Oo = —KoV7B? (123) 


s0; (VFB?) +0: [V7 (VB? vB] + Eo. y76) + veri 


V7 8:6 = -yano . (124) 


Q 


We see that apart from the gradient ¢-term, we obtain another term that involves 
the time-derivative of (/y¢). Hence we need to plug in Equation (43). We rewrite 
the term /0;(,/7¢)/a simplifying the lengthy expression 


fou fou ; A 
GO 79) = ee [V7 (aB — ¢6°)] 
; Ø i 1 il gk pi i 
— Vyko — VIEGAS — aT VION B Okta + — VIB Oi (125) 
=~, [VIBB] + VFB + v7 a8 — yino8! (126) 
Substituting this into (124) yields the modified Faraday’s law (44). 


Appendix C: Derivation of cell centred formulas for FCT 
In the 3+1 decomposition, for the case of a stationary spacetime the induction 


equation can be written in component form as 
O,,/YB? + æl- E.) = 0. (127) 


Integrating this on each of the surfaces bounding a cell with vertexes at 2x} io 


©? J, 24), With | = +1/2, and using the Stokes theorem, we obtain the evolution 


equations for the magnetic flux in CT, for instance 


A® 5 4.4/2,5,k 
a ee = Gi41/2,541/2,k~ Fi41/2,j-1/2,k ~ Fi41/2,5,b4+1/2 + Gi41/2,5,k—-1/2 9 (128) 
where 
Payza = | 2B da? de®, (129) 
= ƏV (zl, ja) 


with each G representing a line integral of the form 


3 
Tk+1/2 


Gi+1/2,j+1/2,k = — f, E3| dz? . (130) 


1 2 
z Zi+1/2°fj+1/2 
k—1/2 
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The fact that each of these integrals appear in the evolution equation of two mag- 
netic fluxes guarantees the conservation of divergence, as will be explained in the 
next Section. 

On the other hand, the numerical fluxes corresponding to the magnetic field com- 
ponents that are returned by the Riemann solver are surface integrals of the electric 


field, for example, the flux in the x?-direction for Bt is 


2592 ia sitiy2 3 i 
AS*F Perce = f, f, Ezs|j41/2 dx” dæ. (131) 


i-1/2 “ Ve-1/2 


The innermost integral is the same as that of Eq. (130), so the average flux can be 
interpreted as 

"aN a! od Ng = AiG; j+1/2,k ; (132) 
where Ĝi j+1/2,k is the mean value of the integral from Eq. (130). To second-order 
accuracy, this integral takes the value Qi j+1 /2,k at the middle of the cell, there- 
fore Gi+1/2,j+1/2,k can be found by interpolating the averaged fluxes from the four 
adjacent cell faces as 


2 p2 2 T2 
A o1 AS*F Lae ASF ETT 
i+1/2,j+1/2,k = J Ar + ATi+1 
E 7 (133) 
A pasa AF yagi 
Ay; Ayj+1 


Since we implemented a cell-centred version of FCT, we are interested in the evolu- 
tion of the average magnetic field at the cell centres. To second order accuracy, the 
rate of change of the average value of the x'—component of the magnetic field is 


d® 
Ti+1/2 dt 


Te+1/2 dỌ dB” Ax; (d® 
—dz = AV;; N 7 
aa dt a a (r 


- : (134) 


Now we substitute Eq. (133) into Eq. (128) and Eq. (128) into Eq. (134). After 
some algebra, we finally obtain eqs. 45 and 46. 


Appendix D: Discretisation of V - B and zero-divergence initial 
conditions 
CT schemes aim to maintain to zero at machine precision the discretisation of the 


divergence given by 


1 


. B ijk = >< 
(V ) „j,k AV; j,k 


(Burnsa — ®i—1/2,j,k + Bi j+1/2,k— 
(135) 


®; j—1/2,k + Bi j,k+1/2 — jaan) ) 
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which can be thought of as the volume average of the quantity 0,(7!/?B*) in the 
given cell. 

When calculating the evolution equation for (V - B); j,k, each of the integrals G 
appear with opposite signs in the expression for d®/dt (128) and cancel to machine 
precision. Therefore, if this discretisation of the divergence was originally zero, it 
will be zero to machine precision during the rest of the simulation. 

However, in the cell-centred version of FCT employed here, we lack information 
concerning the magnetic flux at cell faces, so Equation (135) cannot be used to 
monitor the creation of divergence. We will therefore find a derived quantity that 
we can monitor based on the other available quantities. 


We calculate the average value of the divergence of eight cells sharing a vertex as 


1 
(V. B)i+1/2,j+1/2,k+1/2 = Ay 5 AV(V. B) list, j+l2,k+ls . (136) 
l1,l2,l3=0,1 


When substituting Eq. (135), the right hand side of Eq. (136) consists of a sum of 
terms of the form 


` (Èit3/2,j+l2,k+ls — Paap es + Hp phate — Ďi—1/2,j+l2,k+l3) > 
l2,l3=0,1 


for each direction. Using the same second-order approximation as for the time- 
update, 


A Ax; 
AVi je BE 5 x Zy (Pitt /2,5,4 +4794 4) s (137) 


this becomes 


B!AV 
1+lı 
D [hn 


lı,l2,l3=0,1 ae 


Finally, summing over the three directions, we recover Eq. (47). Since the same 
second-order approximation is used both for the definition and for the time update 
of B°, the definition of divergence given by Equation (47) is conserved to machine 
precision during each evolution step. 

To obtain a divergence-free initial condition, we calculate the magnetic field as 
the curl of a vector potential. First, we calculate the magnetic flux at each of the 
cell faces as 


Pipi /aje = Ait1/2,j+1/2,k — Ait1/2,j—1/2,k ~Ait1/2,j,b41/2 + Ai4s/2,j,h—-1/2 » (138) 
where A are line integrals of the vector potential along the cell edges 


3 
Tk+1/2 


Ai+1/2,j+1/2,k = l. A3| dz? . (139) 


1 2 
Fe Pi41/2°5 41/2 
k—1/2 
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Then we use again the second order approximation from Equation (137) to find 
the average magnetic field components at the cell center. By construction, in this 
way we obtain a divergence-free initial condition using either of the discretization 
of divergence in Eqs. (135) or (47). 
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